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A  procedure  is  developed  to  calculate  the  radiation 
field  as  a  function  of  position,  direction,  and  wavelength 
within  a  spherically  symmetric  circumstellar  dust  shell. 
The  dust  shell  is  assumed  to  consist  of  grey  isotropically 
scattering  dust  particles  in  thermal  equilibrium  with  the 
radiation  field  and  to  be  characterized  by  seven  parameters- 
CD  the  radius  of  the  central  star,  (2)  the  inner  radius  of 
the  shell,  (3)  the  outer  radius  of  the  shell,  (4)  the  total 
optical  depth  of  the  shell,  (5)  an  index  which  specifies 
the  density  distribution,  (6)  the  albedo  of  the  dust  parti- 
cles, and  (7)  the  temperature  of  the  central  star. 

The  temperature  distribution  and  the  radiation  field 
within  several  model  dust  shells  are  determined,  and  used  to 
calculate  for  each  shell  the  spectral-energy  distribution  of 
the  radiation  emitted.   These  models  are  compared  to  those 
of  other  authors.   It  is  found  that  fitting  the  visual 
spectral-energy  distribution  of  an  observed  object  does  not 
necessarily  mean  that  the  model  parameters  so  obtained  are 


consistent  with  the  infrared  part  of  the  objects' spectral 
energy  distribution,  and  that  the  shape  of  the  spectral 
energy  distribution  can  be  strongly  dependent  on  the  index 
in  the  assumed  power  law  density  distribution. 

The  procedure  is  applied  to  the  two  infrared  objects-- 
HD  45677  and  R  Monocerotis. 


SECTION  I 
INTRODUCTION 

This  paper  will  consider  the  problem  of  radiative 
transfer  in  a  spherical  shell  of  dust  particles  illximinated 
by  a  central  star.   In  such  a  shell  a  dust  particle  will 
absorb,  reemit  and  scatter  the  radiation  from  both  the  cen- 
tral star  and  the  other  particles  of  the  shell.   Since  the 
star  will,  in  general,  be  much  hotter  than  the  particles, 
the  radiation  field  will  consist  of  a  visual  part,  due  pri- 
marily to  the  star  and  an  infrared  part  due  primarily  to 
the  dust.   Therefore,  the  solution  of  the  transfer  problem 
must  take  into  account  absorption,  reeraission  and  scattering 
of  both  these  parts  of  the  radiation  field. 

The  solution  of  this  transfer  problem  will  yield  the 
quantities  which  characterize  the  radiation  field  at  any 
point  in  the  shell  and  at  any  frequency  -  the  intensity  and 
the  flux  for  example.   Using  these  quantities  we  can  de- 
scribe the  flow  of  energy  from  the  star  to  the  outer  edge 
of  the  shell  and  predict  the  observational  characteristics 
of  the  shell.   Also,  from  a  theoretical  point  of  view  the 
radiation  field  within  the  dust  shell  is  of  considerable 
interest,  since  it  determines  a  number  of  physical  quanti- 
ties related  to  the  matter  in  the  shell.   The  temperature 
of  the  dust  particles,  the  radiation  pressure  on  the  dust 


and  gas  and  the  excitation  and  ionization  states  of  the  gas 
are  all  dependent  on  the  radiation  field. 

From  a  knowledge  of  the  dust  temperature  within  the 
shell  we  can  limit  the  number  of  possible  materials  of  which 
the  dust  particles  can  be  composed.   According  to  Gilman 
(1969),  at  temperatures  above  1500°K,  refractory  silicates 
such  as  Al2Si05  will  be  the  most  important  condensates 
around  oxygen-rich  stars  while  carbon  will  dominate  around 
carbon-rich  stars  -  below  this  temperature  a  number  of  other 
compounds  can  begin  to  condense.   A  knowledge  of  the  radia- 
tion pressure  as  a  function  of  distance  from  the  star  will 
allow  us  to  compare  its  effect  on  the  material  there  with 
that  of  other  forces  such  as  the  gravitational  force  and  the 
force  due  to  high  velocity  particles  which  might  be  ejected 
from  the  central  star.   The  ionization  and  excitation  state 
of  the  gas  associated  with  a  dust  shell  can  be  of  consider- 
able importance  in  calculating  models  for  collapsing,  rotat- 
ing photostars  interacting  with  an  external  magnetic  field 
(Prentice  and  ter  Haar,  19  71) . 

The  observational  characteristics  of  the  shell  which 
the ^solution  of  the  transfer  problem  would  yield  are  of 
considerable  interest  because  of  the  recent  availability  of 
a  great  deal  of  observational  information  at  infrared  wave- 
lengths.  In  the  past  few  years  a  large  number  of  objects 
have  been  discovered  which  emit  more  infrared  radiation 
than  expected.   Many  of  these  infrared  objects  emit  most  of 
their  radiation  in  the  region  to  the  red  of  ly  and  have 


gone  undetected  in  the  past  because  of  the  low  sensitivity 
of  detectors  in  that  region.   Since  1965  the  region  between 
1  and  20iJ  has  been  surveyed  extensively  and  several  thousand 
infrared  sources  have  been  discovered.   A  partial  list  of 
these  surveys  follows:  the  California  Institute  of  Technology 
survey  at  2.2ii  which  discovered  approximately  5600  sources 
(Neugebauer  and  Leighton,  19  69) ,  surveys  for  red  stars  using 
IN  plates  by  Hetzler  (1937) ,  Chavira  (1967),  Ackermann  et  al. 
(1968) ,  and  Ackermann  (19  70) ,  and  surveys  at  wavelengths 
longer  than  lOy  (e.g.,  Low  1970). 

Many  of  the  objects  discovered  in  these  surveys  have 
an  excess  of  radiation  within  a  given  infrared  wavelength 
interval,  as  compared  with  a  black-body  with  an  effective 
temperature  appropriate  to  the  objects  observed  spectral 
type.   This  excess  of  infrared  radiation  is  usually  referred 
to  as  an  "infrared  excess."   Most  of  the  5600  objects  dis- 
covered in  the  2.2y  survey  have  been  identified  with  late- 
type  giant  stars  (Grasdalen  and  Gaustad,   19  71)  but  about  50 
of  the  reddest  of  these  objects  could  not  be  identified 
with  cataloged  visual  sources.   These  very  red  sources  are 
manifestations  of  a  wide  variety  of  physical  objects,  in- 
cluding according  to  Neugebauer  et  al.  (19  71) ,  HII  regions, 
planetary  nebulae,  M  supergiaht  stars,  late-type  giant 
stars,  RV  Tauri  stars,  novae.  Be  stars,  and  T  Tauri  stars 
and  related  objects.   Most  of  the  infrared  excesses  asso- 
ciated with  the  T  Tauri  stars  and  related  objects  are  as- 
sumed to  result  from  thermal  radiation  from  a  circumstellar 


dust  shell.   This  assumption  is  also  made  for  Be  stars  but 
in  their  case  it  has  been  argued  by  some  authors,  e.g. 
Woolf  et  al.  (19  70) ,  that  the  observed  infrared  emission  is 
consistent  with  free-free  emission  from  ionized  hydrogen. 
For  those  objects  other  than  the  Be  stars  and  the  T  Tauri 
stars  and  related  objects  thermal  radiation  by  dust  is  prob- 
ably present  but  may  not  be  the  principal  cause  of  the 
observed  infrared  excess. 

Since  a  general  solution  of  the  transfer  problem  for 
circumstellar  dust  shells  would  be  readily  applicable  to  the 
T  Tauri  stars  and  related  objects  a  more  detailed  descrip- 
tion of  these  objects  will  be  given  at  this  point. 

The  T  Tauri  stars  and  related  objects  form  an  especial- 
ly interesting  group  from  both  an  observational  and  a 
theoretical  point  of  view.   Poveda  (1965a)  pointed  out  in 
1965  that  T  Tauri  stars,  which  for  the  last  20  years  have 
been  assumed  to  be  young  gravitationally  contracting  objects, 
might  have  thick  circumstellar  dust  shells.   He  argued  that 
such  a  shell  would  be  a  remnant  of  the  contracting  cloud 
from  which  the  star  formed  and  that  it  would  be  bright  in 
the  infrared.   Shortly  after  this  Mendoza  (19  66)  found  that 
most  T  Tauri  stars  do  have  infrared  excesses  and  later 
Mendoza  (1968)  found  that  the  bolometric  luminosity  of  these 
objects  was  from  1.3  to  6.6  times  that  expected  from  the 
visual  observations.   For  R  Monocerotis,  an  object  closely 
related  to  the  T  Tauri  stars,  Mendoza  (19  68)  found  the 
bolometric  luminosity  to  be  58  times  that  expected  from  the 


visual  data.   Figure  9  shows  the  spectral-energy  distribu- 
tion of  R  Monocerotis  and  it  can  be  seen  that  most  of  the 
observed  luminosity  is  emitted  beyond  ly. 

The  T  Tauri  stars  are  Ha  emission  stars,  and  in  the 
early  1940 's  all  such  emission  stars  were  studied  spectro- 
scopically  because  of  the  intrinsic  value  of  their  spectra. 
It  was  not  until  1945  that  Joy  (1945)  recognized  T  Tauri 
stars  as  a  distinct  class  of  emission-line  variables  asso- 
ciated with  nebulosity. 

From  a  study  of  eleven  stars,  Joy  established  the  cri- 
terion to  be  used  in  assigning  stars  to  this  new  class.   He 
describes  in  detail  each  of  the  eleven  stars  and  says  that 
the  spectra  in  general  are  F5-G5,  but  that  there  is  a  small 
variation  in  spectral  type  with  the  brightness.   Also, 
absorption  lines  are  usually  lacking  but  many  bright  lines 
of  low  excitation  are  present.   He  points  out  the  similarity 
between  the  spectra  of  these  stars  and  the  solar  chromo- 
sphere.  For  five  of  the  stars  he  finds  color  excesses,  which 
he  attributes  to  strong  selective  absorption  of  the  nebulosi- 
ty.  He  finds  that  usually  the  emission  lines  seen  are  dis- 
placed to  the  violet  of  the  absorption  lines  (if  any  absorp- 
tion lines  are  present) ,  and  that  it  is  very  difficult  to 
say  anymore  about  the  radial  velocities  of  these  stars. 

In  a  review  article  on  T  Tauri  stars  and  related 
objects,  Herbig  (1960)  points  out  that  the  only  conclusive 
test  for  membership  in  the  T  Tauri  class  is  provided  by  the 
spectroscopic  data. 


The  classification  system  in  use  today  is  described  in 
and  employed  by  the  Third  Edition  of  The  General  Catalogue 
of  Variable  Stars  (Kukarkin  et  al.,  19  69).   The  criteria 
given  in  the  catalogue  are  those  adopted  by  the  lAU  Commis- 
sion 27.   These  criteria  are  spectroscopic  with  one  exception, 
which  is,  that  if  the  variable  is  connected  with  diffuse 
nebulae  it  is  designated  InT  whereas  if  it  is  not  it  is 
designated  IT.   The  prototype  is  stated  to  be  T  Tauri. 

The  catalogue  breaks  the  irregular  nebular  variables 
into  several  classes,  one  of  which,  InT,  represents  the  T 
Tauri  stars.   Its  class  Is  is  the  same  as  Hof fmeister ' s  RW 
Aur  class,  and  differs  little  from  the  InT  class.   Three 
more  of  the  classes  defined  by  the  catalogue,  In,  Ins,  and 
Inb,  contain  young  objects  of  spectral  class  F-M  with  ir- 
regular light  variations  (these  have  been  called  the  Orion 
variables) .   Of  its  nebular  variable  class  the  catalogue 
class  Ina  (early  spectral  type  Orion  variables)  differs  the 
most  from  class  InT. 

In  some  of  the  current  literature  stars  falling  into 
all  these  classes  are  referred  to  as  T  Tauri  stars  and 
related  objects,  although  usually  only  in  the  sense  that 
the  star  is  in  an  early  evolutionary  stage  can  it  be  grouped 
with  the  T  Tauri  variables. 

R  Mon,  the  unusual  object  observed  by  Mendoza,  is  not 
included  in  class  InT  because  it  is  assigned  by  a  Russian 
author  to  spectral  class  A-Fep,  which  is  too  early  a 
spectral  class  for  an  InT  variable.   Joy  (1945)  originally 


classified  it  as  approximately  G,  and  Mendoza  (1968)  classi- 
fies it  as  K.   Thus  this  star  has  been  assigned  spectral 
classifications  from  A  through  K.   This  points  out  an  impor- 
tant fact;  namely,  the  spectral  classification  of  T  Tauri 
stars  is  very  difficult  because  of  the  lack  of  a  strong 
absorption  spectrum. 

Most  of  these  stars  are  rather  faint;  the  brightest  of 
the  class  is  about  m^  =  +  9  at  maximum  light.   The  magnitude 
usually  quoted  for  T  Tauri  variables  is  the  mean  of  the 
observed  magnitude,  even  though  this  may  not  be  the  most 
meaningful  value. 

For  the  last  twenty  years  T  Tauri  stars  have  been 
assumed  to  be  pre-main-sequence  (PMS)  stars  and  recently  PMS 
models  have  become  available  which  might  help  to  verify  this 
assumption.   The  best  proto-star  models  at  present  are  those 
of  Hayashi  (1966)  and  Larson  (19  69b) .  Both  of  these  models 
are  directly  applicable  to  the  problem  of  the  T  Tauri 
phenomena . 

For  his  model  Larson  assumes  that  the  initial  proto- 
stellar  cloud  is  spherical,  not  rotating,  and  free  of 
magnetic  fields.   He  does  not  assume,  as  Hayashi  does,  that 
the  collapse  is  homologous  and  that  the  density  follows  a 
polytropic  law.   Initially  the  cloud  has  the  Jeans  radius 
(R  =  VglGp]"-*-/^,  where  Vg  =  the  velocity  of  sound,  G  =  the 
gravitational  constant,  p  =  the  density) ,  and  a  density  of 
1.10  X  10  ^^   qm/cmr' .      The  time  scale  for  collapse  in  this 
model  goes  as  l//p,  so  that  the  density  gradient  increases 


with  time,  resulting  in  a  dense  core,  which  finally  gains  the 
size  and  mass  of  a  star.   The  core  then  evolves  toward  the 
Main  Sequence  (MS)  while  accreting  matter  from  the  lower  den- 
sity cloud  about  it. 

If  the  mass  of  the  cloud  is  2.5Mg,_or  less,  the  proto- 
star  will  become  visible  somewhere  along  the  lower  part  of 
the  fully  convective  Hayashi  track;  if  it  is  somewhat  more 
massive  it  becomes  visible  at  some  point  on  the  radiative 
track  running  horizontally  toward  the  MS;  if  it  is  very  mas- 
sive (as  massive  as  an  0  or  B  MS  star)  it  does  not  become 
visible  until  reaching  the  MS. 

The  time  required  in  Larson's  model  for  the  accretion 
of  the  envelope  depends  on  the  mass  of  the  proto-star,  and 
ranges  from  10^  to  10   years. 

Hayashi' s  computations  start  with  a  proto-star  of  much 
higher  density,  and  in  his  model  the  accretion  of  the  cloud, 
which  brings  the  proto-star  to  the  beginning  of  the  fully 
convective  phase,  takes  place  in  a  few  decades. 

The  two  models  differ  in  their  treatment  of  the  shock 
phenomena  at  the  interface  between  the  shell  and  the  dense 
core,  but  the  latter  parts  of  their  evolutionary  tracks  are 
in  agreement.   In  both  models  we  have  a  dense  core  and  a 
circumstellar  shell,  and  Strom  (19  72)  points  out  that  we 
should  expect  some  young  stellar  objects  to  show  evidence  of 
the  remnants  of  their  shells.   These  proto-star  models 
neither  explain  nor  predict  the  variability,  line  emission, 
and  mass  ejection  observed  for  the  T  Tauri  stars. 


There  are  a  nuinber  of  reasons  for  believing  that  the  T 
Tauri  stars  are  PMS  objects:  (1)  their  large  red  and  infrared 
excesses  and  deficiencies;  (2)  their  position  on  the  H-R 
diagram  (above  and  below  the  MS);  (3)  their  association  with 
dense  interstellar  matter;  (4)  their  association  with  0  and 
B  stars;  (5)  their  apparent  rotational  velocities;  (6)  their 
anomalous  Li  abundance;  (7)  their  polarization. 

As  noted  earlier,  Joy  (1945)  found  that  almost  half  of 
his  T  Tauri  stars  had  color  excesses.   Since  then  other 
authors  have  found  large  infrared  excesses  for  some  stars, 
and  infrared  deficiencies  for  others. 

Mendoza  (19  66)  did  photometry  in  the  UBVRIJKLM  bands 

from  .36iJ  to  5y  for  2  6  stars,  including  T  Tauri,  RW  Aur,  R 

Mon  and  Lk  Ha  120.   His  observations  show  red  excesses  for 

most  of  the  stars,  a  deficiency  (Ey_j^  -  0.2)  for  DF  Tau  most 

likely  caused  by  the  Ha  line  in  emission  in  the  V  band,  and 

infrared  excesses  (from  0.9  to  5y)  for  all  the  stars 

(E„  „  =  8.5   for  R  Mon) . 
V-M 

Two  years  later  Mendoza  (1968)  published  another  set 
of  photometric  observations  on  Johnson's  nine  color  system; 
this  time  for  33  stars,  some  of  which  he  had  observed 
previously  (Mendoza,  1966) .   His  results  were  similar  to 
those  of  1966.   He  states  that  the  two  stars  with  the  largest 
infrared  excesses  are  R  Mon  and  R  CrA.   Most  of  the  objects 
observed  by  Mendoza  (1966,  19  68)  were  T  Tauri  stars. 

In  making  the  observations  described  above  Mendoza 
used  the  60"  infrared  telescope  at  the  Catalina  observatory 
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in  Arizona.   He  notes  that  making  observations  of  these 
objects  is  made  difficult  by  their  faintness  and  their 
close  association  with  nebulosity.   He  tried  various  dia- 
phram  sizes  from  36  to  13  seconds  of  arc  in  diameter,  and 
found  that  in  many  cases  a  diaphram  smaller  than  13  seconds 
was  actually  needed. 

From  his  observed  fluxes,  Mendoza  (19  68)  computes 
bolometric  corrections,  and  then  absolute  magnitudes  for 
several  of  these  objects.   He  finds  that  they  have  high  lu- 
minosities, and  speculates  that  Joy  (1945)  classified  them 
as  dwarfs  because  of  their  rather  wide  absorption  lines. 

To  explain  his  19  66  observations  Mendoza  (1966) 
advanced  two  theories:  (1)  a  multiple  star  system  in  which 
one  of  the  components  is  a  very  cool  infrared  star,  and  (2) 
a  core-envelope  model  in  which  the  UBV  photometry  pertains 
to  the  core  and  the  infrared  photometry  to  the  envelope.  In 
this  paper  (Mendoza,  1966)  he  does  not  favor  one  theory  over 
the  other. 

In  the  second  paper  Mendoza  (1968)  states  that  his 
previous  core-envelope  model  is  probably  the  explanation  for 
the  large  infrared  excesses  observed  for  objects  such  as  R 
Mon  and  R  CrA.   For  the  objects  with  small  infrared  excesses 
(infrared  just  a  few  times  the  visual)  he  suggests  that  the 
cause  is  heating  of  the  stellar  wind  by  cosmic  rays. 

He  points  out  that  his  observations  indicate  that  the 
envelope  around  the  core  behaves  like  a  neutral  filter  in 
the  visual  (UBV)  region.   Assuming  that  all  the  infrared 
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radiation  comes  from  this  envelope,  he  finds  that  it  would 

have  a  radius  of  about  10  AU  and  a  mass  of  about  10"-^  M 

® 

Low  and  Smith  (19  66)  pioblished  the  results  of  their 
photometry  of  R  Mon  at  20y.   They  found  that  at  that  wave- 
length the  star  was  about  one  fourth  as  bright  as  Mondoza 
(1966)  found  it  to  be  at  3.4y. 

To  explain  their  observations  they  computed  the 
spectral-energy-distribution  to  be  expected  from  an  optical- 
ly thin  circiimstellar  dust  shell  in  which  scattering  and 
reabsorption  of  infrared  radiation  emitted  by  the  dust  can 
be  neglected.   The  theoretical  distribution  they  obtained 
matched  the  available  observations  fairly  well. 

Strom  et  al.  (1972)  observed  42  stars  in  NGC  2264  at 
1.6,  2.2,  and  3.4y,and  found  infrared  excesses  and  some 
infrared  deficiencies.   The  authors  state  that  these  obser- 
vations are  an  indication  that  the  stars  have  circumstellar 
shells.   They  used  profiles  of  the  observed  hydrogen  lines 
to  estimate  the  surface  gravities  of  some  of  the  PMS  A  stars 
in  their  sample  of  42  NGC  2264  stars.   From  these  surface 
gravities  they  computed  luminosities,  and  found  that  those 
stars  that  seemed  the  faintest,  as  compared  with  their 
computed  luminosities,  all  had  infrared  excesses.   They  sug- 
gest that  these  infrared  excesses  are  due  to  a  dust  shell, 
radiating  in  the  infrared,  energy  it  has  absorbed  from  the 
star  in  the  visual  part  of  the  spectrum.   The  stars  that 
they  found  to  have  the  largest  infrared  excesses,   were  those 
that  lie  below  the  MS. 
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Strom  et  al.  (19  72)  also  found  that  these  shells  have 
a  very  high  ratio  of  total  to  selective  absorption,  which 
would  suggest  large  dust  particles.   The  dust  shells  of 
lowest  temperature  were  around  the  latest  stars. 

They  speculate  that  the  infrared  deficiencies  observed 
could  be  due  to  a  very  cool  dust  shell  at  a  great  distance 
from  the  central  star.   Of  course,  observations  in  the  far 
infrared  should  find  the  missing  radiation. 

Martin  Cohen  (1973b) has  observed  stars  in  NGC  2264,  IC 
5146,  VI  Cyg  and  NGC  7000  (The  North  American  Nebulae).   NGC 
7000  contains  the  T  Tauri  star  Lk  Ha  190,  which  was  discover- 
ed by  Herbig  in  1958,  and  has  recently  brightened  visually 
by  6  magnitudes. 

Cohen  used  the  60"  Catalina  infrared  telescope,  and 
made  observations  in  all  wavelengths  in  a  single  night  (this 
is  very  important  since  these  objects  vary  rapidly) .    He 
concluded  that  for  most  of  these  stars  (many  of  them  not  T 
Tauri  stars)  thermal  emission  from  circumstellar  dust  is  the 
cause  of  the  infrared  excesses  he  observed.   In  a  few  cases, 
he  says  that  free-free  emission  from  a  circumstellar  shell 
may  be  important  out  to  about  lly  . 

Breger  (1972) ,  like  Cohen,  finds  evidence  for  circum- 
stellar shells  about  PMS  stars  in  NGC  2264.    As  well  as 
infrared  excesses,  he  finds  light  variability,  which  he  at- 
tributes to  circumstellar  shells. 

Strom  et  al.  (19  71)  found  a  small  amount  of  continuous 
and  Ha  emission,  similar  to  that  found  for  T  Tauri  stars,  in 
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the  spectra  of  some  ft-F  stars  in  young  clusters.   They  attri- 
bute this  to  very  thin  circumstellar  shells.   They  also  say 
that  the  strength  of  these  features  increases  as  the 
circumstellar  absorption  does. 

According  to  Strom  (19  72)  we  should  expect  to  see  free- 
free  (electron-H  ion)  and  free-bound  emission  from  envelopes 
in  the  infrared,  if  they  have  sufficient  size  and  density. 
Dust  envelopes,  as  stated  previously,  would  show  thermal 
infrared  reradiation;  thus,  Strom  (19  72)  ascribes  the  cir- 
cumstellar emission  in  the  optical  region  of  T  Tauri  spectra 
to  gas  (this  is  the  emission  that  washes  out  the  normal 
absorption  spectrum) .   He  says  that,  in  a  joint  paper,  Kuhi 
and  the  Stroms  indicate  that  gas  shell  emission  may  also  do- 
minate the  infrared  region. 

In  support  of  his  idea  that  gas  and  not  dust  shell 
emission  is  dominant  in  the  infrared,  Strom  (19  72)  offers 
the  following  argument  and  observational  evidence.   He 
finds  that  for  some  young  stars  (especially  A-F  stars)  with 
circumstellar  shells.  Ha  emission  and  infrared  excesses 
there  exists  a  correlation  between  log  (Ha  intensity)  and 
the  L  magnitude  (3.4y).   In  an  optically  thin  case  of  the 
dust  model  the  dust  emission  is  proportional  to  the  optical 
depth  (t  =  p  x  length) ,  whereas  the  gas  emission  (Ha) 
depends  on  p^  x  length.   Therefore,  unless  we  believe  that 
all  circumstellar  shells  have  the  same  densities,  the  dust 
theory  is  difficult  to  defend,  because  we  would  not  get  the 
observed  correlation. 
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According  to  Herbig,  in  the  Taurus  cloud,  Orion 
nebulae,  NGC  2264,  NGC  6530,  and  IC  5146,  we  find  that  if  we 
fit  the  brightest  T  Tauri  stars  to  a  MS  with  a  distance  modu- 
lus of  6  the  Ml  to  M3  stars  observed  are  too  bright.   Thus, 
the  observed  MS  in  these  regions  appears  to  be  turned  up  on 
the  red  end;  however,  the  observed  MS  does  not  curve  upward 
more  and  more  as  we  observe  later  and  later  stars,  but  runs 
parallel  to  the  Zero  Age  Main  Sequence  (ZAMS)  after  sharply 
breaking  off  from  it  at  about  AO .   This  observed  MS  forms  a 
band  with  its  upper  edge  several  magnitudes  above  the  ZAMS. 

Herbig  (1960)  believes  that  the  uncertainties  in  the 
observations  and  theory  are  large  enough  so  that  we  need  not 
be  concerned  by  the  fact  that  the  observed  MS  does  not  curve 
upward.   He  also  believes  that  part  of  the  spread  in  the 
observed  lower  MS  is  caused  by  a  spread  in  the  times  of 
stellar  formation  in  the  cluster. 

Hayashi's  evolutionary  tracks  for  proto-stars  would 
lead  us  to  expect  young  stars  to  lie  above  the  ZAMS. 

Poveda  (196  5a)  first  brought  attention  to  objects 
falling  below  the  ZAMS  in  NGC  2264.   He  said  that  their  lo- 
cation below  the  ZAMS  was  caused  by  thick  circumstellar 
dust  shells. 

Strom  et  al.  (1972),  as  we  have  seen,  were  able  to 
show  that  those  stars,  in  NGC  2  264,  lying  below  the  MS  would 
be  above  it,  if  their  circiimstellar  shells  were  removed. 

Presumably,  the  stars  below  the  MS  are  younger  than 
those  above  it,  and  have  not  had  time  to  flatten  their 
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circumstellar  shells.   According  to  Mendoza  (196  8) ,  the  T 
Tauri  stars  with  the  least  infrared  excess  are  the  oldest, 
and  have  the  most  flattened  disk-like  shells.   If  we  are 
looking  at  a  disk  edge  on,  we  should  see  a  large  discrepan- 
cy between  the  luminosity  as  determined  from  the  spectral 
classification,  and  the  bolometric  luminosity  observed. 
Strom  (19  72)  finds,  for  a  sample  of  10  Herbig  Ae  and  Be 

stars,  that  for  three,  L„^  >  >  L,  ^t^„^.^-^  . 

'   sp      DOlometric 

According  to  Strom  et  al.  (19  72) ,  the  presence  of 
late-type  stars  near  the  MS  in  young  clusters  is  caused  by 
circumstellar  obscuration,  and  not  by  a  large  spread 
(several  x  10^  years)  in  times  of  stellar  formation,  as 
argued  by  Herbig  (1960) . 

T  Tauri  stars  are  found  directly  associated  with  the 
clouds  in  which  they  presumably  formed.   They  can  not  be 
interlopers  from  the  general  field  around  the  region  in 
which  they  are  found,  since  their  space  density  averages  1 
or  2  orders  of  magnitude  above  that  of  the  field  stars  of 
the  same  luminosity  (Herbig,  19  60)  . 

Cohen  (1973c)  has  observed  the  T  Tauri  type  objects 
at  the  tips  of  several  cometary  nebulae.   Mendoza  (19  66) 
states  that  T  Tauri  and  R  Mon  illuminate  bright  nebulosities. 
Both  of  these  nebulosities  are  variable  cometary  nebulae. 
R  Mon  is  in  NGC  22  61,  Hubbies  variable  nebulae,  and  T  Tauri 
is  in  NGC  1555,  Hind's  nebulae. 

Poveda  (196  5c)  describes  the  cometary  nebulae 
(elephant  trunks,  bright  rims,  etc.)  as  HII  regions  expand- 
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ing  into  a  surrounding  HI  region,  or  as  a  HII  region  expand- 
ing past  an  area  of  high  density,  compressing  and  deforming 
it.   If  either  of  these  explanations  is  correct,  the  tip  of 
such  a  nebulae  will  be  a  region  of  high  density,  suitable 
for  star  formation.   Thus,  cometary  nebulae  are  very  young 
phenomenon  illuminated  by  even  younger  stars. 

If  we  assume  the  stars  have  the  space  velocities  of 
the  eddy  of  gases  from  which  they  formed  (about  2  km/sec) , 
and  if  we  estimate  the  size  in  km  of  the  tip  of  the  nebulae, 
we  can  estimate  how  long  it  would  take  the  star  to  drift 
out  of  the  tip,  and  thereby  arrive  at  an  upper  estimate  for 
the  age  of  the  star.   Poveda  does  this  for  several  nebulae, 
and  finds  ages  running  from  5  x  10-^  to  50  x  10   years. 

Mendoza  (19  68)  also  makes  a  calculation  to  obtain  the 
ages  of  some  T  Tauri  stars.   He  uses  the  PMS  evolutionary 
tracks  of  Iben  to  determine  the  age  of  a  star  from  its 
position  on  the  H-R  diagram.   The  ages  Mendoza  obtains  in 
this  way  agree  with  those  of  Poveda  for  some  stars  and 
disagree  for  others. 

For  R  Mon  the  two  methods  agree  rather  well,  but 
for  the  other  stars  Poveda 's  results  are,  in  general,  lower 
than  Mendoza' s. 

Herbig  (I960)  states  that  in  regions  such  as  Orion,  we 
find  T  Tauri  stars  associated  with  0  and  B  stars,  which  are 
considered  to  be  young  objects.   It  seems  logical  then,  to 
consider  the  T  Tauri  stars  to  be  of  similar  age  (hence 
young),  but  of  lower  mass,  and  thus  not  as  greatly  evolved. 
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Joy  (1945)  classified  the  T  Tauri  stars  as  MS  dwarfs, 
which  we  now  know  to  be  incorrect.   Those  T  Tauri  stars 
that  do  show  absorption  lines,  have  very  wide  lines,  which 
might  indicate  that  they  were  luminosity  class  V  objects. 
However,  it  is  now  believed  that  these  abnormally  wide 
absorption  lines  are  caused  by  rotational  broadening  (Herbig, 
1960). 

The  rotational  velocities  found  for  these  stars 
(20  km. sec  <  v  sin  i  <  65  km/sec)  are  unsually  high  for 
stars  of  their  spectral  class.   Normal  G-K  stars  have 
V  sin  i  <  15  km/sec. 

If  we  compute  the  radii  of  the  T  Tauri  stars  from 
their  observed  luminosity  and  color- temperature,  and  then  as- 
sume they  contract  without  loss  of  angular  momentum  along  a 
Hayashi  track,  we  can  find  their  rotational  velocities  on  the 
MS.    When  this  is  done,  the  computed  values  are  in  close 
agreement  with  observed  values  of  v  sin  i  for  MS  stars  of 
comparable  mass.    It  may  of  course  be  true  that  the  line 
broadening  is  not  caused  by  rotation. 

The  T  Tauri  stars  as  a  group  appear  to  be  overabundant 
in  Li,  as  compared  with  the  Sun,  by  1  or  2  orders  of  magni- 
tude, which  is  taken  as  an  indication  of  their  youth  (Herbig, 
1960) .   Herbig  compared  the  Li  abundance  of  the  T  Tauri 
stars  to  that  of  chondritic  meteorites,  which  presumably 
have  the  Li  abundance  of  the  pre-solar  nebulae. 

It  is  possible  that  the  large  abundance  of  Li  is  being 
caused  by  Li  production  in  the  photospheric  layers  of  the  T 
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Tauri  stars,  in  which  case  it  may  not  be  an  indication  of 
their  youth.   It  seems  unlikely  though  that  such  Li  produc- 
tion is  actually  occuring,  since  the  most  plausible 
mechanisms  suggested  for  Li  production  require  large,  strong 
magnetic  fields,  which  have  not  been  observed. 

Zellner  (19  70) ,  Capps  and  Dyck  (19  72)  and  Breger  and 
Dyck  (19  72)  observe  polarization  for  some  stars  in  the  near 
and  far  infrared.   Zellner 's  measurements  for  R  Mon  indicate 
polarization  typical  of  large  grains  (radius  "^   O.ly)  while 
the  other  authors  cited  indicate  that  their  measurements 
favor  electron  scattering  similar  to  that  observed  for  Be 
stars.   In  a  more  recent  paper  Zellner  and  Serkowski  (19  72) 
state  that  R  CrA  exhibits  polarization  similar  to  R  Mon. 
They  suggest  that  in  objects  such  as  R  Mon  and  R  CrA  the 
dust  cloud  producing  the  infrared  excess  is  in  the  form  of 
an  equatorial  disk  and  the  light  seen  at  shorter  wavelengths 
is  scattered  light  which  has  filtered  out  the  poles  of  this 
disk. 

Breger  and  Dyck  find  that  between  3000  and  7500  A, 
four  out  of  a  sample  of  35  stars  in  NGC  2264  show  intrinsic 
polarization.   Most  of  the  stars  they  observed  were  not  T 
Tauri  stars,  but  of  the  four  showing  polarization  one  is  a 
T  Tauri  star.   Breger  and  Dyck  comment  that  its  polarization 
has  a  different  A  dependence  from  that  they  observe  for  the 
other  three  polarized  stars  in  their  sample,  and  that  its 
polarization  is  also  different  from  that  of  T  Tauri,  as 
given  by  Serkowski  (19  71) . 
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These  investigations  indicated  the  presence  of  circum- 
stellar  shells  about  T  Tauri  stars  and  stars  in  clusters  as- 
sociated with  T  Tauri  stars,  such  as  NGC  2264. 

Some  other  important  characteristics  exhibited  by 
young  stars  and  T  Tauri  stars  are  the  following:  (1)  mass 
loss  and  the  P  Cyg  effect,  (2)  mass  accretion  and  the  inverse 
P  Cyg  effect,  (3)  ultraviolet  excess,  and  (4)  correlations 
of  the  visual  magnitude  and  the  infrared  emission,  of  the  Ha 
UV  and  IR  emissions,  and  between  mass  infall  and  the  V 
magnitude. 

Most  T  Tauri  stars  show  a  P  Cyg  type  spectrum,  that  is, 
their  emission  lines  have  blueshif ted  absorption  features 
superimposed.   According  to  Strom  (19  72)  and  Herbig  (19  60) , 
the  theory  that  this  effect  is  caused  by  an  expanding  envel- 
op with  self  absorption  is  widely  believed. 

Herbig  (1960)  states  that  the  star  T  Tauri  has  a  set 
of  two  absorption  lines  superimposed  on  its  H  and  Ca  emis- 
sion lines,  corresponding  to  radial  velocities  of  -70  and 
-170  km/sec,  and  that  RY  Tau  has  one  line  corresponding  to  a 
radial  velocity  of  -90  km/sec. 

From  an  analysis  of  the  forbidden-line  radiation,  com- 
ing from  the  small  forbidden-line  region  surrounding  T  Tauri, 
Herbig  says  that  Varsavsky,  in  1960,  finds  a  mass  loss  of 
about  .5  X  10"^  Mgj/yr  for  T  Tauri,  assuming  that  the 
forbidden-line  region  is  750  AU  in  diameter,  and  that  the 
region  is  sustained  by  material  rising  at  170  km/sec. 

The  accretion  of  matter  by  some  young  stars  is  indicated 
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by  the  presence  of  an  inverse-  P  Cyg  spectrum  (sometimes 
called  a  YY  Orionis  spectrum) .   Out  of  a  sample  of  2  3  UV- 
excess  stars  observed  by  Walker  (1969)  in  NGC  2264,  he  finds 
10  show  an  inverse  P  Cyg  effect.    None  of  these  are  T  Tauri 
stars,  but  they  are  all  closely  related  to  T  Tauri  stars.  In 
19  72,  Walker  (19  72)  reports  finding  two  more  young  stars 
with  inverse  P  Cyg  effect,  RR  Tauri  and  e^C  Orionis. 

T  Tauri  spectra  show  a  strong  continuum  in  the  visual, 
which  is  sometimes  called  the  'blue  continuum',  and  a 
continuum  starting  in  the  AX3700-3800  region  and  rapidly 
rising  toward  the  shorter  wavelengths.   Unlike  the  blue 
continuum,  this  UV  continuum  adds  an  appreciable  amount  of 
energy  to  the  luminosity  of  the  star.    According  to  Herbig 
(1960) ,  this  UV  continuum  has  been  attributed  to  synchrotron 
emission  like  that  from  some  solar  grains,  and  to  the  running 
together  of  the  higher  Balmer  lines.   The  latter  of  these 
hypotheses  he  refers  to  as  Bohm's  hypothesis.   The  higher 
Balmer  lines  blend  together  because  of  turbulent  broadening, 
according  to  this  hypothesis. 

Mendoza  (1966)  finds  UV  excesses  for  the  majority  of 
the  26  stars  he  observes.   Anderson  and  Kuhi  (1969)  give 
detailed  observations  of  the  T  Tauri  star  AS  209  (1900:  16^^ 
43"*. 6,  -14°  13'),  and  report  that  it  shows  a  large  UV  excess. 
They  say  they  favor  Bohm's  hypothesis,  since  Kuhi's 
observed  line  width  for  the  H3  line,  indicates  a  turbulent 
velocity  as  large  as  100  km/sec,  which  would  be  more  than 
enough  to  blend  the  higher  Balmer  lines. 
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In  general.  Ha  emission  and  UV  excess  in  T  Tauri  stars 
is  correlated,  but  for  individual  stars  there  does  not  appear 
to  be  a  correlation  between  the  two  as  the  star  varies  in 
brightness.   In  the  case  of  AS  209,  Anderson  and  Kuhi  looked 
for  such  a  correlation,  and  found  none. 

Various  correlations  have  been  noted  between  observed 
T  Tauri  phenomena.   Herbig  finds  that  as  RW  Aur  gets  fainter 
in  the  V  band,  it  gets  redder.   Anderson  and  Kuhi  (1969)  find 
that  AS  209  increases  its  UV  excess  as  it  gets  fainter  in  the 
V  band.   Cohen  (19  73a)  finds  that  as  T  Tauri  gets  fainter  in 
the  visual  it  gets  brighter  in  the  near  infrared. 

Walker  (19  69)  does  a  detailed  study  of  SU  Ori,  an  Inb 
star  with  inverse  P  Cyg  spectrum  and  UV  excess,  and  finds 
that  when  the  redward  displaced  absorption  feature  is  most 
prominent,  the  star  is  brightest.   He  shows  five  spectro- 
grams of  SU  Ori  taken  when  the  star  was  at  different  magni- 
tude, and  the  redward  displaced  absorption  line  disappears 
as  the  star  gets  dimmer. 

Poveda  (1965b,  1965c)  suggests  that  the  circumstellar 
shells  about  T  Tauri  stars  and  related  objects  are  pre- 
planetary  systems.   He  says  that  some  of  the  variability  of 
these  stars  must  be  due  to  eclipses  by  proto-planets  (large 
clouds  along  the  orbits  of  the  planets-to-be) ,  and  that 
phenomena  such  as  FU  Orionis  occur  because  the  circumstellar 
material  can  form  into  large  particles  very  rapidly  (a  few 
years) ;  thus,  rapidly  dropping  the  optical  depth  of  the 
material,  and  brightening  the  star. 
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In  light  of  the  observed  characteristics  of  the  T 
Tauri  stars  and  related  objects  described  above  it  is  not 
surprising  that  there  have  been  a  number  of  papers  published 
in  the  past  few  years  dealing  with  the  problem  of  radiative 
transfer  in  circumstellar  envelopes.   In  each  of  these 
papers  the  author  makes  a  number  of  assiomptions  about  the 
object  or  objects  he  is  attempting  to  model  in  order  to  sim- 
plify the  radiative  transfer  problem.   These  simplifying  as- 
sumptions tend  to  reduce  the  accuracy  and  the  versatility  of 
the  resulting  model  and  should  therefore  be  avoided  if 
possible. 

In  all  of  these  papers  it  has  been  assumed  that  the 
envelope  is  spherically  symmetric  about  the  central  star  and 
that  the  effects  of  the  spherical  symmetry  are  significant. 
If  the  entire  shell  were  metrically  thin  and  distant  from 
the  central  star  it  would  be  possible  to  simplify  the  prob- 
lem by  assuming  the  shell  to  be  stratified  into  parallel 
planes.   This  simplification  could  also  be  made  if  the  outer 
layers  of  the  envelope  were  optically  thick.   However,  it 
appears  from  the  observations  that  the  circumstellar  clouds 
about  the  T  Tauri  stars  and  related  objects  can  not  be 
adequately  represented  by  a  model  based  on  a  plane  parallel 
approximation.   The  assumption  of  spherical  symmetry  while 
much  better  than  the  simple  plane  parallel  assumption  is 
still  an  idealization  since  it  has  been  pointed  out  that 
visually  these  objects  appear  to  have  a  very  complicated 
outer  structure  (e.g.,  Herbig,  1968). 
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In  all  these  papers  it  is  assumed  that  a  power  law  of 
the  form  K(r)  =   CC/r'^]  governs  either  the  voliime  absorption 
coefficient  or  the  number  density  of  particles.   Here  r  is 
the  distance  from  the  central  star,  C  is  a  constant,  and  n 
is  a  variable  index  (Chandrasekhar ,  19  34)  .   Some  assumption 
as  to  the  functional  relationship  between  K  and  r  is  neces- 
sary (see  Section  II) ,  and  this  one  seems  to  allow  reason- 
able flexibility  with  the  adjustment  of  only  one  parameter 
(the  constant  is  determined  by  the  geometry  and  total 
optical  depth  of  the  shell) .   In  a  few  of  these  papers  it  is 
further  assumed  that  n  =  0;  that  is,  that  K(r)  is  constant 
with  r  . 

Several  authors  (Stein,  1966;  Low  and  Smith,  1966; 
Krishna  Swamy,  19  70)  have  considered  the  less  general  prob- 
lem of  radiative  transfer  in  an  optically  thin  circumstellar 
shell  in  which  scattering  and  reabsorption  of  infrared  radia- 
tion emitted  by  the  dust  can  be  neglected.   Stein's  paper 
pertains  to  the  observed  infrared  excesses  of  certain  B 
stars;  Low  and  Smith's,  to  the  much  larger  infrared  excess 
of  R  Monocerotis  and  Krishna  Swamy 's,  to  the  infrared  excesses 
of  several  infrared  objects,  among  them,  VY  Canis  Majoris. 

Stein  in  his  paper  states  that  the  total  optical  depth 
of  the  shells  he  is  considering  is  about  10   .   For  a  shell 
this  optically  thin  it  is  quite  reasonable  to  assume  that 
the  star  is  the  only  significant  energy  source  in  the  system. 
In  the  paper  by  Swamy  and  in  the  the  one  by  Low  and  Smith, 
the  objects  modeled  (e.g.,  VY  CMa  and  R  Mon)  show  an  amount 
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of  infrared  radiation  as  compared  with  the  optical  which 
indicates  that  these  objects  have  optically  thick  shells 
(Mendoza,  1966;  Low  et  al.,  1970). 

Huang  in  one  paper  (Huang,  1969a)  considers  the  two 
extreme  cases  of  very  small  and  very  large  optical  thickness 
and  then  in  a  second  paper  (Huang,  19  69b)  considers  the  more 
general  problem  of  intermediate  optical  thickness.   In  both 
these  papers  Huang  makes  the  assiamption  that  the  opacity  of 
the  dust  is  proportional  to  the  geometrical  cross  section 
of  the  grains. 

In  the  first  paper  he  solves  the  transfer  problem  for 
an  optically  thin  shell  in  much  the  same  manner  as  Stein 
(1966) .   His  approach  to  the  problem  of  the  optically  thick 
shell  is  analogous  to  that  of  Chandrasekhar  (19  34)  for 
extended  stellar  atmospheres. 

In  the  second  paper  he  solves  the  transfer  equation 
for  a  spherical  shell  having  its  inner  boundary  in  contact 
with  the  surface  of  the  central  star,  assuming  grey 
particles,  the  Eddington  approximation,  isotropic  scattering 
and  Local  Thermodynamic  Equilibrium.   He  does  not  actually 
state  that  the  shell  is  in  contact  with  the  surface  of  the 
star  but  that  the  radiation  at  the  inner  boundary  of  the 
shell  is  completely  diffuse.   For  all  practical  purposes  as 
long  as  the  inner  edge  of  the  shell  is  within  a  few  stellar 
radii  of  the  surface  of  the  star  the  radiation  field  can  be 
considered  diffuse. 

Huang  (19  71)  published  a  third  paper  dealing  specif i- 
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cally  with  shells  distant  enough  from  their  central  stars 
such  that  the  stars  could  be  considered  point  sources.   He 
points  out  in  this  paper  that  the  assumption  of  a  complete- 
ly diffuse  radiation  field  throughout  the  shell  represents 
an  idealized  case  of  one  extreme  while  the  assumption  that 
the  star  is  a  point  source  represents  the  other  extreme. 

In  his  second  paper  Huang  (19  69b)  argues  that  the  as- 
sumption that  the  particles  scatter  radiation  as  grey-bodies 
is  valid.   He  states  that  in  the  case  of  interstellar  red- 
dening blue  light  is  preferentially  scattered  out  of  the 
line  of  sight,  whereas  in  the  case  of  a  circumstellar  dust 
cloud  we  observe  scattered  radiation  coming  in  all  direc- 
tions; such  that,  the  amount  of  radiation  of  a  certain  wave- 
length which  is  lost  in  direct  transmission  along  the  line 
of  sight  is  regained  due  to  the  overall  scattering  process. 

For  his  intermediate  case  Huang  (19  69b)  assumes  that 
the  radiation  field  can  be  broken  into  two  parts  -  a  visual 
part  and  an  infrared  part.   He  then  proceeds  to  solve  the 
overall  transfer  problem  as  though  it  were  two  less  diffi- 
cult problems:  (1)  a  nonconservative  scattering  problem  in- 
volving only  the  stellar  radiation  in  the  visual  (the 
infrared  emission  of  the  star  is  assumed  to  be  negligible) 
and  (2)  the  problem  of  an  atmosphere  in  local  thermo- 
dynamical  equilibrium  for  the  infrared  radiation  emitted  by 
the  dust,  with  an  energy  source,  namely,  the  visual  radia- 
tion absorbed  in  the  nonconservative  optical  region.   There- 
fore, the  results  for  his  intermediate  case  described  in  his 
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second  and  third  papers  (Huang,  1969b,  1971)  refer  to  these 
two  broad  wavelength  regions . 

In  his  results  Huang  (19  71)  states  that  the  radiation 
field  and  the  temperature  distribution  do  not  critically 
depend  upon  the  index  n  in  the  expression  assumed  for  the 
density  distribution  and  indirectly  that  the  ratio  of  the 
inner  to  the  outer  radii  of  the  shell  is  an  unimportant  fac- 
tor.  In  Section  IV  it  is  shown  that  these  results  are  not 
in  general  correct. 

Huang  does  not  apply  this  method  to  any  observed 
object,  but  in  1970  Herbig  (1970)  uses  it  successfully  in 
the  development  of  a  model  for  VY  Canis  Majoris,  and  recent- 
ly Apruzese  (19  74)  has  used  it  to  model  the  infrared  object 
HD  45677  (a  Be  star  with  a  large  infrared  excess) .    The 
approach  followed  by  Apruzese  parallels  Huang's  for  the  in- 
termediate case  and  diffuse  radiation  except  that  he 
modifies  Huang's  technique  slightly  to  allow  the  determina- 
tion of  the  radiation  field  in  ten  spectral  regions  to  the 
blue  side  of  ly  .   For  each  of  these  spectral  regions  he 
has  in  effect  a  problem  of  nonconservative  scattering. 
Apruzese 's  model  will  be  considered  further  in  Section  IV. 

Larson  (1969b)  in  studying  the  emission  of  collapsing 
protostars  finds  an  approximate  solution  for  the  radiative 
transfer  in  a  spherical  dust  shell  in  which  the  particles 
scatter  none  of  the  radiation  but  absorb  an  amount  propor- 
tional to  X~P,  where  X   is  the  wavelength  and  p  is  a  varia- 
ble index.   He  derives  an  approximate  expression  for  the 
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temperature  distribution  in  the  cloud  and  using  this  solves 
for  the  emitted  spectrum.    The  weighting  function  integra- 
tion method  he  develops  to  compute  the  emitted  spectrum 
from  the  known  temperature  distribution  is  also  used  by 
Herbig  in  the  paper  mentioned  above.   For  the  index  n  which 
is  used  to  establish  the  density  distribution  in  the  cloud 
Larson  uses  1.5  and  states  that  this  is  the  value  determined 
by  his  dynamical  calculations  and  published  in  a  previous 
paper  (Larson,  1969a) . 

Hummer  and  Rybicki  (19  71)  have  a  procedure  involving 
iteration  on  the  Eddington  factor,  K/J,  which  they  use  to 
solve  the  transfer  problem,  under  the  assumption  of  the  con- 
servative grey  case,  for  a  spherical  cloud  extending  from  a 
point  of  radiation  at  its  center  out  to  some  finite  radius, 
R. 

None  of  the  papers  mentioned  above  adequately  deal 
with  the  problem  of  radiative  transfer  in  circumstellar  en- 
velopes of  intermediate  optical  depth.   The  only  two  proce- 
dures that  come  close  to  solving  this  case  are  those  due  to 
Huang  (1969b,  1971),  Larson  (1969b),  and  Apruzese  (1974). 
They  all  assume  that  the  density  in  the  shell  is  proportion- 
al to  l/r'^  and  Huang  and  Apruzese  assume  that  the  particles 
are  isotropic  scatters  -  Larson  neglects  scattering.    The 
first  of  these  assumptions  seems  as  good  as  any  other  for 
the  density  distribution  and  is  not  critical  to  the  solution 
of  the  problem.   The  second  assumption  is  necessary  in 
order  to  take  advantage  of  the  spherical  symmetry  of  the 
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shell.   It  allows  us  to  reduce  the  problem  to  one  involving 
only  two  spacial  coordinates.   If  the  particles  are  assumed 
to  scatter  according  to  an  arbitrary  phase  function  the 
problem  becomes  much  more  complex.   Huang  assiomes  that  the 
opacity  of  the  cloud  is  due  to  the  geometrical  cross  section 
of  the  particles  and  therefore  that  it  is  grey,  whereas, 
Apruzese  and  Larson  assume  that  the  opacity  is  wavelength 
dependent.   The  procedure  described  by  Huang  is  restricted 
to  cases  where  the  central  star  can  be  considered  a  point 
source  or  the  radiation  is  completely  diffuse.   It  is 
further  restricted  to  cases  where  the  Eddington  approxima- 
tion is  valid.   Most  importantly  his  procedure  does  not  de- 
scribe the  radiation  field  at  a  particular  wavelength  but 
merely  breaks  it  into  two  regions  -  visual  and  infrared.  The 
procedure  described  by  Apruzese  is  based  on  Huang's  proce- 
dure and  involves  one  major  improvement  -  his  procedure  can 
be  used  to  determine  the  radiation  field  at  a  number  of 
wavelengths  in  the  visual.   Larson's  procedure,  as  well  as 
neglecting  scattering,  uses  a  rather  crude  approximation 
for  the  temperature  distribution  in  the  shell.   His  proce- 
dure does  however  allow  the  determination  of  the  spectral- 
energy  distribution  emergent  from  the  shell  in  the  visual 
and  infrared. 

In  Sections  II  and  III  a  technique  is  developed 
capable  of  solving  the  transfer  problem  for  the  radiation 
field  at  any  point  in  a  circumstellar  shell,  for  any  wave- 
length, and  any  optical  depth. 
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This  technique  is  applicable  to  cases  where  the 
central  star  can  be  considered  neither  a  point  source  nor  a 
diffuse  source  of  radiation.    To  simplify  the  problem  the 
following  assumptions  are  made:  (1)  the  cloud  is  spherically 
symmetric  about  the  central  star,  (2)  the  cloud  is  made  of 
small  particles  which  have  no  preferential  orientation,  (3) 
the  particles  reemit  absorbed  radiation  isotropically  and 
scatter  radiation  isotropically,  (4)  the  particles  have  an 
albedo  which  is  independent  of  wavelength,  (5)  the  particles 
radiate  as  grey-bodies  at  a  temperature  T  where  T  is  a  func- 
tion of  the  radiation  field  about  the  particle,  (6)  the 
volume  absorption  coefficient  as  a  function  of  the  distance 
from  the  central  star,  K(r) ,  is  constant  with  wavelength  and 
proportional  to  the  density  p(r),  and  (7)  the  optical  depth 
T  and  the  distance  from  the  central  star  are  related  by  the 
usual  expression  (Chandrasekhar,  19  34) , 

dt   ^   r 
dT  =  -  b  ;^  ;  t  =  ^  , 

where  b  is  a  constant. 

This  technique  characterizes  each  dust  shell  by  a  set 
of  parameters  related  to  the  size  of  the  shell  and  the 
density  distribution  within  the  shell,  the  size  and  tempera- 
ture of  the  central  star,  and  the  albedo  of  the  dust 
particles  within  the  shell.   The  equation  of  transfer  is 
solved  in  integral  form  by  a  computer  code  using  an  itera- 
tive procedure. 


SECTION  II 
BASIC  RELATIONS 

In  this  section  the  physical  and  mathematical  relations 
necessary  to  solve  the  transfer  equation  for  the  radiation 
field  in  a  circumstellar  dust  cloud  will  be  developed.   These 
relations  will  be  in  terms  of  a  number  of  physical  parameters 
that  characterize  the  dust  cloud. 

The  parameters  to  be  used  can  be  placed  into  three 
groups,  one  that  specifies  the  geometry  of  the  cloud,  one  that 
specifies  the  albedo  of  the  particles,  and  one  that  specifies 
the  temperature  of  the  central  star.   The  first  group  consists 
of  five  parameters,  the  radius  of  the  central  star,  the  inner 
radius  of  the  cloud,  the  outer  radius  of  the  cloud,  the  total 
optical  depth  of  the  cloud,  and  an  index  which  specifies  the 
density  distribution.   These  five  quantities  will  be  desig- 
nated respectively  as  R^,  Rj ,  R,  Tj ,  and  n.   The  second  and 
third  groups  have  only  one  member  each:   the  albedo  of  the 
particles,  designated  w^,  in  the  case  of  the  second  group  and 
the  temperature  of  the  central  star,  designated  T*,  in  the 
case  of  the  third. 

For  some  set  of  these  seven  parameters  the  dust  cloud 
model  should  be  capable  of  reproducing  the  observed  spectral- 
energy  distribution  of  an  object  such  that  that  set  of  param- 
eters can  be  considered  to  characterize  the  observed  object. 
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Consistent  with  the  simplifying  assumptions  set  forth 
in  Section  I,  relations  can  be  derived  for  a  number  of  useful 
quantities.   Using  the  expression  for  dt  given  in  Section  I 
we  can  derive  the  relation  for  the  absorption  per  unit  volume 
at  a  distance  r  from  the  central  star,  K(r).   We  have 

dr  =  Rdt 

bR^~^dr 

dx  =  -  ^ — ' 


but 


and  therefore. 


dx    = 

K(r)dr, 

K{r) 

bR^-1 

or 


K(r)  =  A 

r 


(1) 


where  B  is  a  constant.  We  can  also  derive  expressions  for 
the  optical  depth  at  r,  T(r) ,  and  for  B  in  terms  of  R^,  R, 
T-r,  and  n.   In  the  following  n  ?^  1  unless  specified  otherwise, 


I 
We  have 


r 

X  (r)  =  /  K(r)dr 

o 


r 
X (r)  =  B/  r""dr 
o 


and,  therefore, 

^.1-n 


(r)  =  B(  ^p:^  )  +  C   ,  (2) 


where  C  is  a  constant  of  integration.   Since  x(Rj)  -  x^  and 
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T (R)  =  0  we  have 

_l-n 
1-n 


TT  =  B  (  5^—  )  +  C  (3) 


„l-n 

0  =  B  ( )  +  C  (4) 

1-n 


from  which  we  obtain, 

C  =  -  B  (  5__  ). 
1-n 

Substituting  this  into  (3)  we  have  for  t^. 


R  1-n         1-n 

.   =  B(  -^ )  -  B(  5 ) 

I       1-n  1-n 


^  ,^  1-n    T,l-n, 
I    1-n 


T^  =  —  (Rj"-"  -  R-  ") 


So 


Tj(l-n) 


B  = 


1-n  _  pl-ri' 
'I 


and  we  obtain  from  equation  (1) 

K(r)  =      ^I^^""^      ,  if  Rj  <  r  <  R     (5) 
(r1-^  -  R  l-^)r" 


=  0,  if  r  >  R,  or  r  <  Rj. 

and  from  equation  (2) 

,„l-n    r-^~^\ 

T(r)  =  -^ •  (6) 

(R^-^  -  R  1-") 
I 
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In  the  case  with  constant  density  we  have  n  =  1  and 
the  relations  for  K(r)  and  T(r)  become 

K(r)  =  T-|-[£n(R/Rj)r]   ,  Rj  <  r  <  R 

=  0,  r  >  R,  r  >  Rj  (7a) 

T(r)  =  T^Jln(R/r)  [  n(R/Ri)]""'"  (7b) 

Equations  (6)  and  (7b)  apply  only  for  R,  <  r  <  R.   For 
r  <  Rj,  T(r)  =  T  ,  and  for  r  >  R,  T(r)  =  0. 

The  spectral-energy  distribution  is  observed  in  terms 
of  flux  received  at  the  surface  of  the  earth  per  unit  wave- 
length, so  this  is  the  quantity  we  want  to  calculate  with 
our  model.   If  we  define  F(A,r)  as  the  flux  at  wavelenth  X 
through  a  unit  area  of  the  shell  of  radius  r  about  the 
central  star,  where  F(X,r)  is  in  units  of  ergs/cm'^  -   k   - 
str  -  sec,  we  have 

F(X,D)   =   F(X,R)R^/D^  , 

where  D  is  the  distance  between  the  object  and  the  earth. 
Note  that  if  D  is  not  a  known  quantity  we  can  still  find  a 
scaled  flux  and  thus  determine  at  least  the  shape  of  the 
spectral-energy  distribution. 
We  also  have 

F(A,r)  =  /  I(A,r,u))cose  du  ,  (8) 

0) 

where  I  (A,  r,  to)  is  the  specific  intensity  of  radiation  at 
wavelength  A,  at  any  point  a  distance  r  from  the  central 
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star,  and  coming  from  some  direction  w  radians  from  the  cen- 
tral star.   The  integration  is  over  47t  steradians.    Because 
of  our  assumption  of  spherical  symmetry  we  can  define  I 
uniquely  with  reference  to  only  two  of  the  three  spherical 
coordinates,  r  and  G;  there  will  be  no  (j)  dependence.   It  is 
also  true  that  I(A,  r,  6)  =  I(A,  r,  -9),  so  the  integration 
over  6  need  only  be  carried  out  from  0  to  it.    Thus, 

F(X,r)  =  27r  /  I(A,r,e)sine  cosG   de 
o 

where  9  is  an  angle  in  a  single  plane. 

In  order  to  determine  F(X,r)  and  the  F(A,D)  we  must  be 

able  to  calculate  I(X,r,6)  numerically.   This  can  be 

expressed  as 

I(A,r,9)  =  /  "•  e-t  S  (A,rg  ,9^)  dt,  (9) 

o 

where  S(A,r  , 9  )  is  the  source  function  for  radiation  of 
wavelength  A  ,  emitted  in  a  direction  9g  radians  from  the 
central  star,  at  any  point  a  distance  r^  from  the  central 
star.   The  integration  is  carried  out  over  optical  depth 
along  the  line  of  sight,  from  a  point  at  a  distance  r  from 
the  central  star  and  in  a  direction  9  radians  from  the  cen- 
tral star.   The  upper  limit  for  the  integration,  t-^, 
represents  the  maximum  optical  depth  attainable  along  the 
line  of  sight  and  will  always  be  less  than  2t-j-.   The  quan- 
tities rg  and  9   are  functions  of  r,  9  and  t.   Here  t  is 
optical  depth  along  the  line  of  sight  and  the  variable  of 
integration.   Thus  we  have  r  =  rg(9,r,t)  and  9^  =  9  (9,r,t), 
The  expression  for  I(A ,r, 0) , \given i in  equation  (9), 
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does  not  help  us  unless  we  can  evaluate  S{X,rQ,Q^) .      The 
source  function  can  be  expressed  as  the  ratio  of  the  emission 
to  the  absorption  coefficient,  that  is, 

S(A,rg,e^)  =  j(A,r^,ee)/K{rg) , 

where  j(A,re,ee)  gives  the  emission  per  unit  volume  at  a 
distance  r^  from  the  central  star,  at  wavelength  X,    and  in  a 
direction  9   radians  from  the  central  star.   Since  the  dust 
model  is  based  on  the  idea  that  the  dust  scatters  part  of  the 
radiation  it  intercepts  and  absorbs  and  reradiates  the  rest, 
we  can  write  S{A,re,6g)  in  the  following  form: 

S(A,r3,63)  =  [JB(A,r^,ee)  +  Jg^^'^e'^e)]  /  ^^^e^     '  ^^^^ 

where  k(X,rg,ee)  has  been  broken  into  two  parts,  scattered 
emission  and  thermal  or  grey-body  emission,  designated 
respectively  jg  and  jg  . 

It  is  evident  that  both  emission  terms  in  equation  (10) 
are  related  to  the  radiation  field,  I(X,re,e),  at  r^.  There- 
fore, we  need  to  determine  the  exact  functional  relationships 
expressing  JB(A,re,ee)  and  Jg(X,rg,eg)  in  terms  of  I(A,re,e). 
In  doing  this  we  will  first  need  to  derive  expressions  for 
^e  "  rg(9,r,t)  and  9^  =  9^(9, r,t). 

From  consideration  of  the  geometry  in  Figure  1,  we 

have 

1/2 
r  {e,r,t)  =  (r^  +  s^  -  2rs  cos9)  ^  ,        (11) 
e 

where  s  is  the  metric  distance  along  the  line  of  sight  to  the 
point  of  emission  and 
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0   =  sin"-*-  [  (r/r„)sine] . 

The  variable  of  integration  along  the  line  of  sight,  t, 

can  be  evaluated  from  the  expression 

s 
t  =  /  K(r)  ds  (12) 

o 

once  we  have  assigned  a  value  to  n  in  the  formula  for  K(r) . 

The  amount  of  radiation  scattered  by  a  particle  will 

be  that  amount  which  is  intercepted  by  the  cross  sectional 

area  of  the  particle  and  not  absorbed;  that  is/  the  amount 

intercepted  times  co  .   The  amount  of  this  radiation  scattered 

in  a  direction  6^  will,  in  general,  be  given  by 

Js(^/rg,6g)  =  Ckpa)y4TT)    /   P  (a)  I  ( A,rg,  9)  d(o,  {13a) 

where  k   is  the  cross  sectional  area  of  the  particle  and 
P(a)  is  a  phase  function  defined  such  that 

(1/47T)/  P(a)d(o  =  1  . 

The  quantity  a  is  the  angle  through  which  the  radiation  is 

scattered  and  is  thus  a  function  of  6  at  r^  and  of  6   (see 

e         e 

Figure  1) . 

Since  we  are  assuming  isotropic  scattering,  where  P(a) 
=  1  for  all  a,  we  can  drop  6   from  our  notation  and  rewrite 
equation  (13a)  as 

Js(A,r3)    =    (kpajQ/47T)    /    I(X,re,e)    dco    .  (13b) 

IX) 

Equation  (13b)  gives  the  radiation  scattered  per  particle 
and  we  want  the  radiation  scattered  per  cm^.   To  derive  the 
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needed  form  we  must  replace  k   in  equation  (13b)  with  the 

ir 

cross  sectional  area  of  all  the  particles  in  a  cm-^,  which  we 
will  denote  k.   Then, 


JsC^.rg)  =  (kWo/4TT)  /  I(A,rg,0)da) 
which  because  of  spherical  symmetry  we  can  rewrite  as 

Js{A,rg)    =    {kcoQ/2)     /    I  (A, r^,  9)  sine    dG  .       (14) 
o 

Two  additional  equations  for  jg  are  formed  by  separat- 
ing I(A,rg,9)  into  a  term  due  to  the  cloud  and  a  term  due  to 
the  star.   Let  y{Q)    be  defined  such  that  y{Q)    =  1  if 

0  <  e  <  e*(rg)  and  yO)  =  0  if  6  >  e*(rg),  where  e*(rg)  = 

_1  R* 
sin   (— )  is  the  angular  radius  of  the  central  star  as  seen 

e 
from  a  distance  r  ;  then 

I(A,r3e)  =  Ic(A,r3,e)  +  Y(e)  I*(A)  . 

In  the  above  equation  I  ,  the  term  due  to  the  cloud,  is  given 

c 

by 

Tc   -t 
I^,(A,r3,e)  =  /   e  ^  S(A,r)  dt 

o 
where  r  =  r(6,rg,t)  is  obtained  analogously  to  r^  (see  equa- 
tion 11)  , 

^c  =  ^m  -  ^(^)   (^b  "  V 
and  Tjj  represents  the  optical  depth  along  the  line  of  sight 
to  the  inner  boundary  of  the  cloud.   The  term  due  to  the  star 
is  I^(A)  =  e~'^bB(A,T*)  ,  where  B(A,  T*)  gives  the  Planckian 
emission  of  the  star. 
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From  equation  (14)  and  the  above  expression  for  I  we 
have 

IT 

Js(^'re)  =  (kajQ/2)  [  /  I^  (A,re,  9)  sin6  dG  + 
o 

/  I^(A)sine  de  ] ,         (15a) 
o 

and  if  we  are  only  interested  in  the  bolometric  emission,  we 

have 


jg(rg)  =(ka)Q/2)[  /  I^(re,e)sine  dS  + 


c 
o 


9*   Tg-Tj   ^m   4 

/   e      (— =^)sin9  d9  ],  (15b) 


TT 


where 


Io(^^'9)  =  /  l(A,r  ,9)dA,and  a  (15c) 

c  c  e 

o 

is  the  Stefan-Boltzmann  constant. 

When  9^  is  very  small  we  can  use  the  following 

approximate  versions  of  equations  (15a)  and  (15b) : 

TT 

Js^^'^e^  =  (^(^0/2)  /  I^(A,r^,9)sin9  d9 
o 

+  (kcjQ/4TT)exp(Te-Tj)R2TrB(X,T^)/r^  (16a) 

and 


jg(rg)  =  (kcoQ/2)  /  I^(rg,9)sin9  d9 
o 

+  (ka)Q/4Tr)exp(Tg-Tj)R2  aT^/r^  (16b) 

where  t   is  the  optical  depth  from  the  surface  of  the  cloud 
to  rg  and  the  other  quantities  used  have  been  previously 
defined. 
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We  now  turn  to  the  problem  of  determining  the  grey- 
body  emission  of  the  particles,  jg(X,rg,9).   As  long  as  any 
elemental  volume  is  small  enough  for  every  particle  in  it 
to  experience  nearly  the  same  radiation  field,  the  amount  of 
energy  absorbed  will  be  independent  of  the  shape  of  the 
element  and  its  orientation.   This  absorbed  energy  will  be  a 
function  of  I  and  the  effective  cross  sectional  area  of  the 
particles  per  unit  volume,  k.   Thus,  we  can  derive  an  expres- 
sion for  k  by  considering  the  energy  absorbed  by  an  idealized 
unit  volume;  namely  a  small  cylinder  of  unit  height  and  unit 
cross  sectional  area  oriented  such  that  the  incident  radia- 
tion is  perpendicular  to  one  end. 

For  such  a  cylinder  the  amount  of  the  intensity  lost 

_x 
due  to  absorption,  I-j^,  will  be  I_  =  1(1  -  e   ),  where  x  here 

is  the  total  optical  depth  along  the  axis  of  the  cylinder, 

or  if  we  assume  the  medium  to  be  thin  enough  so  that  occulta- 

tion  of  one  particle  by  another  are  rare,  we  have 

•^L  "^    ^^part^^cy^  ■'■'  ^^^^^  %art  ^"'^  ^cy  respectively  denote 
the  cross  sectional  areas  of  the  particles  and  the  cylinder. 
Now  we  can  write 

A   ^  =  (1  -  e~'^)A^„  , 
part  '  cy  ' 

and  since  x  <  <  1,  we  have 

A     =  [1  -  (1  -  T)]A   =  AX  , 
part  '^^  ^y 


A   ^  =  A^^K(r)  h  , 
part    cy  ^  ' 

where  h  is  a  length  over  which  we  have  an  optical  depth  x  and 
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K(r)  has  been  assiamed  constant  over  this  length.   For  a  cyl- 
inder of  unit  volume  we  can  take  A_...  =  1  and  h  =  1;  then 
A  ^j,^/unit  vol  =  K(r)  ,  and  k(r)  =  K(r). 

We  now  need  to  derive  a  relation  between  the  effective 
cross  sectional  area  of  the  particles  in  a  volume  and  the 
effective  surface  area  enclosing  that  volume.   Let  A  be  the 
sum  of  the  effective  surface  areas  of  all  the  individual 
particles  in  a  volume  element  V  and  let  k  be  the  sum  of  their 
effective  cross  sectional  areas.  Define  the  effective  cross 
sectional  area  of  a  single  particle,  to  be  its  geometrical 
cross  sectional  area  averaged  over  all  possible  orientations. 
Now  if  the  number  of  particles  in  V  is  large,  and  they  have 
no  preferential  orientation;  then  k  will  not  vary  with  the 
direction  of  observation  o),  and  we  will  be  able  to  derive  a 
relation  for  A  in  terms  of  k;  thus  defining  A. 

For  the  volume,  V,  the  energy  absorbed,  E^,  is  given 
by  the  relation 


E^  =  /  I(co)k  do)  , 

0) 


if  the  density  is  low  enough  such  that  occuUations  of  one 

particle  by  another  can  be  neglected.   The  energy  emitted, 

4 
Eg,  is  given  by  Eg  =  aT  A.   Therefore, 


/  I{a))k  do)  =  aT'^A  , 


0) 


if  the  particles  are  in  thermal  equilibrium  with  the  radia- 
tion field.   This  can  be  rewritten  as 
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T^  =  (k/Aa)  /  Kw)  do)  ,  (17) 

and  if  we  let  W  =  k/A,  we  have 

T^  =  (w/a)  /  Kco)  do)  .  (18) 

(0 

This  equation  is  true  for  any  volume  of  particles  and  any 
radiation  field. 

If  we  apply  equation  (18)  to  a  volume  situated  in  a 
black-body  cavity  at  temperature  Tg,  we  have 

Tg  =  (W/a)  /  B(T3)  dco, 

and  since  B(Tg)  is  constant  with  u),  we  have 

Tg   =    [W47rB(Tg)  ]/a    .  (19) 

Since  B(Tg)  is  not  a  function  of  A,  it  represents  the 
Planckian  integrated  over  A.   Thus, 

CO 

B(Tb)  =  /  B(A,Tg)dA  =  (aTj)/7r  , 
o 

and  equation  (19)  becomes  T**  =  4WTg,  from  which  we  obtain 

W=  1/4,  a  value  indentical  to  that  obtained  for  a  sphere. 

Substituting  this  value  for  W  into  equation  (18)  gives 

us 

1/4 
T  =[(l/4a)  /  I(a))da)] 

This  equation  can  be  used  to  find  the  temperature  of  any 
grey-body  due  to  an  arbitrary  radiation  field,  I((jj), 
provided  that  the  body  is  tumbling  randomly. 

Since  W  is  the  ratio  of  the  effective  cross  sectional 
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area  to  the  effective  surface  area  of  a  particle,  A  =  4k 
for  any  particle. 

From  our  assumptions  that  the  particles  have  no  pref- 
erential orientation  and  that  they  are  in  thermal  equilib- 
rium with  the  radiation  field,  we  have  that  A  =  4k  ,  and 

P 

that 

F(A,T)  =  (1  -  w^)TiB(A,T)  , 

where  A  is  the  effective  surface  area  of  a  particle,  T 
is  the  equilibrium  temperature  of  the  particle,  B(A,T)  is 
the  Planckian  and  F(X,T)  is  the  flux  per  unit  effective 
area  of  the  particle.   The  thermal  isotropic  emission  of  a 
single  particle  can  now  be  written  as 

Jb(^)  =  (A/4Tr)F(A,T)  , 


Jb(^)  =  kp(l-cOo)B(A,T)  . 

The  corresponding  emission  for  a  unit  volume  of  such 
particles  will  be 

JB(A,r^)  =  K(rg)  (l-a)Q)B(X,T)  . 

Note  that  the  temperature  T,  like  K,  is  a  function  of  the 
distance  from  the  central  star,  r. 

If  this  emission  is  considered  only  bolometrically , 


F(T)  =  (l-u)Q)aT^  , 


and 
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j^(r  )  =  [K(r^)/7r]  (l-a)o)aT4  .  (20) 

The  condition  of  thermal  equilibri\im  gives  us  the 
relation 

E^(r)  =  Eg(r)  ,  (21) 

where  E  (r)  and  E  (r)  denote/ respectively, the  energy  absorbed 
and  emitted  by  a  unit  volume  of  particles  at  a  distance  r 
from  the  central  star.   Now  from  equations  (20)  and  (21)  we 
have 

E^(r)  =  /  (1-w^)  (a/Tr)K(r)T'*(r)du  , 

a.  O 

(0 

and  since  the  emission  is  isotropic 

E^(r)  =  4K(r)  (l-u^)  aT^^  (r)  . 

From  this  formula  for  E  we  can  obtain  the  following 
expression  for  the  temperature, 

T(r)  =  {[l/4a(l-a)o)]  [E^(r)/K(r)]}f/^  , 

For  the  non-bolometric  case  we  have, 

E  (r)  =  K(r)  (l-co  )  //  I  (X  ,r ,  6)  dXdw  , 

giving  us 

1/4 
T(r)  =  [(l/4a)  //  I(A,r,e)dAdco] 

Because  of  spherical  symmetry  we  can  write  this  as 

^  1/4  , 

T(r)  =  [(TT/2a)  //  I(A,r,e)sine  dX  d0] 

oX 
If  we  divide  the  radiation  field  into  two  parts,  as 
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we  did  before,  we  have 

T(r)  =  {{l/4a)[2Tr  /  I  (r,e)sine  de  + 
o  ^ 


9*  aT*4   T-T,         1/4 
/   —Z —  e    sine  de]  }    ,     (22) 


where  I^(r,e)  is  the  bolometric  intensity  defined  by  equation 
(15c) .   When  the  central  star  can  be  considered  a  point 
source  we  have  the  additional  equation, 

IT 

T(r)  =  {(l/4a)  [27T  /  I  (r,e)sine  d9  + 
o  "" 

e    ■"  (R^aT-^/r^)]} 

Note  that  in  equation  (9)  I  is  a  functional  of  S  and 
that  in  equation  (10a)  S  is  a  function  of  j .   We  have  seen 
that  j  is  a  functional  of  I;  therefore,  I  is  a  functional  of 
S  and  S  is  a  functional  of  I.   Thus,  equations  (9)  and  (10a) 
can  be  rewritten  as  a  recurrence  formula  for  I(A,r  ,a)),  which 
can  be  used  in  an  iterative  solution  for  the  radiation  field 
throughout  the  cloud  and  hence  for  the  desired  spectral- 
energy  distribution. 

The  recurrence  formula  is 

Ij^(X,r,u)  =  /    e~^  S[  /  I^_-,^(A,re,a))da)]dt,        (23) 

O  0) 

where  the  subscripts  n  and  n+1  refer  to  iterations  n  and 
n+1.   The  true  radiation  field  consistent  with  a  given  set 
of  parameters  will  be  given  by 

I(X,r,aj)  =  lim  :^(A,r,a)) 
n->-<» 
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This  completes  the  definition  and  derivation  of  the 
basic  quantities  needed  to  implement  an  iterative  solution 
of  equation  (11) . 

Once  the  proper  set  of  the  parameters,  t  ,  R  ,  R  and 
n,  is  determined,  we  can  obtain  values  for  the  density  with- 
in the  cloud,  p(r),  and  the  total  mass  of  the  cloud  M. 

3 
If  we  let  D  equal  the  number  of  particles  per  cm  and 

d  equal  the  average  diameter  of  a  particle,  where  d  <  <  R, 

then 

o 

D  =  — CSA/cm =  particles/cm^  . 

CSA/particle 

The  cross  sectional  area  of  the  particles  per  cubic  centi- 
meter is  a  function  of  r  and  is  given  by  K(r) ;  thus, 

D(r)  =  ^^4^   =  I  K(r)d-2  . 
2 

For  the  volume  and  mass  of  a  particle,  V  and  m  respective- 

P      P 

ly,  we  have 


V  =  —  77  (^)   cm  /particle, 
P    3    2 


and 


m  =  p  Y"  (j)   grams/particle  , 
where  p   is  the  specific  gravity  of  a  grain.   The  density  of 

IT 

the  cloud  at  a  distance  r  from  the  central  star  will  then  be 
given  by 


p(r)  =  D(r) 


mp 
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2 


p(r)  =  i  PpK(r)  (|)  . 

Substituting  for  K{r) ,  we  find 

P(r)  =  4  Pn  ^  t^i ^ ]  ^  •  (24) 

3   P    Rl-n_Rl-n  P 

From  the  above  expression  for  p(r)  we  can  derive  an 
expression  for  the  total  mass  of  the  cloud.   We  have 
R 


'I 


M  =  /   (4TTr^)  p^^Mr  , 


from  which  we  obtain, 

8Trp„T,{l-n)  (r3-^^  -  R  3-n)^ 

M  =  P-J^ . i .        (25) 

3(3-n) (Rl-n  -  R  1-n) 


SECTION  III 
NUMERICAL  PROCEDURE 

In  this  section  a  procedure  based  on  the  equations  of 
Section  II  is  developed,  to  effect  a  numericsl  solution  for 
the  radiation  field  in  the  dust  shell.   This  procedure  will 
involve  iterative  numerical  integrations  of  equation  (23) 
to  obtain  a  set  of  convergent  series  of  XS^,    r-,co)  .   Here  i 
specifies  the  i^"  such  set  and  r^^,  therefore  specifies  the 
distance  from  the  central  star  at  which  the  I(A,r-,a))   for 
that  set  are  to  be  evaluated.   The  n  signifies  that  this  is 
a  value  obtained  in  the  n''-'^  iteration.  All  the  integrations 
will  be  performed  using  Gauss-Legendre  formulas. 

In  order  to  accomplish  the  necessary  integrations  it 
will  first  be  necessary  to  establish  an  integration  grid  in 
keeping  with  the  boundary  conditions  imposed  by  the  geometry 
of  the  shell  and  the  central  star.   Consider  the  spherical 
shell  of  dust  extending  from  R-j-  to  R  and  divided  into  N 
thinner  concentric  shells  each  with  the  same  optical  thick- 
ness, T  /N  (see  Figure  2  ) .   The  inner  and  outer  radii  of 
these  shells  are  a  function  of  t  and  can  thus  be  found  from 
the  equation 

r(T)  =  [r1-^  -  (T/Tj)  (Rl-n  -  rj^""")]  ^         (26) 
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For  each  of  the  N  points  lying  on  the  semi -diameter  of 
the  cloud,  SR,  (see  Figure  2)  at  distance  r-  from  the  central 
star,  we  assign  M  values  of  the  angle  9,  which  will  define  M 
lines  radiating  from  each  point.   We  then  break  these 
radial  lines  into  segments  of  approximately  equal  optical 
depth,  where  the  end  points  of  the  segments  are  all  the 
points  at  which  the  radial  lines  intersect  the  N  circles  of 
radii  r^. 

Thus,  about  each  of  the  N  points  on  SR  we  have  M 
radial  lines,  and  each  of  these  lines  is  broken  into  from  1 
to  2N  segments,  where  the  actual  number  of  segments  for  any 
particular  line  depends  on  how  many  of  the  N  circles  of 
radii  r^  that  line  intersects. 

The  M  values  of  9  used  at  each  of  the  N  points  along 
SR  would  in  general  be  the  Gauss-Legendre  points  for  the 
interval  0  to  tt;  however,  it  is  usually  advantageous  to 
break  this  interval  into  two  intervals  with  M/2  Gaussian 
points  each  -  one  from  0  to  9jj^(ri)  and  the  other  from  Q^   to  tt, 
where  9ni{ri)  is  an  angle  that  must  be  determined  for  each  of 
the  N  points  r^.   If  r-  is  small  enough  so  that  the  central 
star  shows  an  appreciable  disk,  then  9jjj(rj^)=9*  (rjL)  ,  but  if  rj^ 
is  large  enough  so  that  the  star  can  be  considered  a  point 
source  9jjj(rj^)  will  be  an  angle  between  0  and  it/2.  The  value 
of  9  (rj^)  in  this  case  should  be  determined,  such  that,  it 
lies  between  Tr/2  and  the  value  of  9  at  which  I(r^,9)  has  its 
greatest  value.   Defining  9^(rj^)  in  this  manner  allows  us  to 
concentrate  our  Gaussian  points  on  the  part  of  the  interval 
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0  to  ir  where  1(9)  is  varying  most  rapidly.   In  one  special 
case  the  M  Gaussian  points  must  be  determined  in  a  different 
manner.   When  ri  =  r^  =  R  the  M  9's  are  restricted  to  the 
interval  0  to  7r/2  since  for  Tr/2  <  0  <  tt,  I(R,e)  obviously 
vanishes. 

We  have  broken  each  of  the  M  radial  lines  into  segments 
as  a  first  step  in  the  determination  of  the  locations  of  L 
Gaussian  points  to  be  used  in  the  numerical  evaluation  of 
equation  (23).   Figure  (2)  shows  a  few  representative  radial 
lines  for  the  i*^  of  the  N  points  on  SR.   Also,  the  radial 
lines  are  shown  broken  into  segments  with  end  points  a 
distance  Sj^  from  SR,  where  s^  is  measured  along  the  radial 
line  and  k  is  an  index  ranging  from  1  to  2N. 

Let  X  and  y  be  the  rectangular  coordinates  of  a  point 
in  a  right-handed  coordinate  system,  where  the  central  star 
is  the  origin,  and  the  semi-diameter  SR  is  the  positve 
y-axis.   Also,  let  the  slope  and  y-intercept  of  a  radial  line 
be  denoted  by  m  and  b  respectively;  then,  each  pair  of 
values  for  m  and  b  defines  a  particular  radial  line.   In 
terms  of  previously  defined  quantities,  m  =  Cot  (6)  and 
b  =  r^^  .   When  b  and  m  are  specified,  the  rectangular  coordi- 
nates of  all  the  points  where  a  particular  radial  line  inter- 
sects the  N  circles  of  radii  r.  can  be  found  by  solving,  for 
X  and  y,  the  system  of  equations 

x^  +  y^  -  r^^  =  0 

y  -  mx  -  b  =  0,  for  i  =  1  to  N. 
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We  are  only  interested  in  points  to  one  side  of  SR,  so  we 
exclude  all  solutions  with  x  >  0.   Points  which  are  occulted 
by  the  star  are  also  excluded. 

The  starting  point  of  a  radial  line  will  lie  on  SR, 
and  have  coordinates  x  =  0  and  y  =  r.;  the  k   point  of 
intersection  on  that  line  will  have  coordinates  x^   and  y-.  . 
Thus,  the  distance  along  the  line  from  the  starting  point  to 
the  k^  point,  which  is  denoted  s,  in  Figure  2  ,  is 

Sk  =  [Xk2  +  (r._y^)2pl/2  , 

To  accomplish  the  integration  of  equation  (2  3)  we  must 
first  perform  the  integration  indicated  by  equation  (12) 
along  each  of  the  radial  lines  using  the  end  points  of  the 
segments  on  these  lines  as  the  limits  of  integration. 
Equation  (11)  is  used  here  to  find  the  distance  from  the 
central  star  to  a  point  on  the  line  and  hence  the  values  of 
K(r)  needed  by  (12)  from  (5)  or  (7a). 

Prom  (12)  we  will  obtain  for  each  of  the  N  x  M  lines  of 
sight  (LOS)  a  set  of  accurate  optical  depths,  tj^,  and  cor- 
responding distances,  Sj^,  measured  along  the  LOS.   Because 
of  the  discontinuity  caused  by  the  central  region  of  the 
dust  cloud  (r  <  RNR)  it  is  helpful  to  break  the  LOS  into  two 
intervals  to  consist  of  L/2  Gaussian  points. 

Unless  the  LOS  intersects  the  central  region  of  the 
cloud  (r  =  r  =  RNR) ,  the  midpoint  of  each  LOS,  at  which  we 
can  break  the  LOS  into  two  segments,  is  roughly  determined 
by  dividing  the  total  number  of  points  by  two.   However,  if 
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tile  shell  r-|^  is  crossed  two  "midpoints"  are  established 
(say,  MIDI  and  MID2)  where  the  first  such  point  is  the  point 
of  intersection  of  the  LOS  and  r,  closest  to  the  starting 
point  and  the  second  is  the  next  point  further  along  the  LOS. 
The  segments  from  the  beginning  of  the  LOS  to  MIDI  and  from 
MID2  to  the  intersection  of  the  LOS  with  r,  are  to  be 
considered  the  two  halves  of  the  LOS  (see  Figure  2:) .   One 
additional  case  remains;  if  the  LOS  crosses  the  shell  r^  and 
then  intersects  the  surface  of  the  star  (r  =  R*)  the  LOS  is 
to  be  considered  terminated  at  r-^    .   This  will  obviously 
OGciar  whenever  9(ri)  <  e*(r-)  for  a  LOS. 

Next  each  half  of  the  LOS  is  broken  into  19  segments 
(this  is  the  specific  number  used  in  the  present  version  of 
the  computer  code  -  it  can  be  varied  if  necessary)  as  nearly 
the  same  length  in  optical  depth  as  computing  time  will 
allow  (the  computer  code  is  given  a  parameter  which  sets  the 
precision  desired) .   This  is  done  by  an  iterative  procedure 
in  which  the  computer  code  uses  an  Aitken-Lagrange  inter- 
polation routine  (Hildebrand,  1956;  Mullen,  1974)  to  inter- 
polate in  the  table  of  s^^  and  tj^  values  for  the  set  of  Sj^ 
values  corresponding  to  a  new  set  of  evenly  spaced  t-^   values. 
A  Gaussian  integration  routine  is  then  used  to  evaluate 
equation  (12), thus  redetermining  the  values  of  tj^  .   If  the 
values  of  t,  found  by  this  integration  are  not  equally 
spaced  (to  the  precision  requested)  another  interpolation 
and  integration  are  performed,  and  this  process  is  continued 
until  the  desired  accuracy  in  spacing  is  achieved. 


52 


If  the  original  set  of  sj^  and  t-^   values  contains  only 

two  points  (the  case  where  the  LOS  starts  at  the  outer 

boundary  r  and  then  intersects  r  without  intersecting 

r  _-, )  the  procedure  described  above  breaks  down  since  a  set 

of  two  points  can  have  no  midpoint.   In  this  special  case 

the  code  sets  up  a  more  extensive  table  containing  20  points 

(midpoint  at  10)  by  breaking  the  interval  S2  to  s-^  into  19 

equal  segments  thus  obtaining  20  Sj^  values  and  then  finds 

the  corresponding  t,  values  by  integration.   The  resulting 
ic 

tables  of  20  s-^   values  (equally  spaced)  and  t-^   values  (un- 
equally spaced  unless  n  =  0  in  density  relation)  are  then 
used  in  the  iterative  process  to  find  the  set  of  t,  values 
at  equal  spacing  in  optical  depth. 

Since  we  want  to  evaluate  equation  (2  3)  by  Guassian 
quadrature  we  must  be  able  to  evaluate  the  integral,  at  the 
Gaussian  points,  t]^„,  on  the  interval  0  to  Xj^j    therefore, 
we  next  interpolate  in  our  table  of  Sj^  and  tj^  values  for 
the  s,  values,  designated  Sj^   corresponding  to  the  t^ 
values. 

We  have  that  S(t)  =  S(T(t)),  where  t  is  the  optical 
depth  measured  from  the  outer  boundary  of  the  cloud  to  the 
point  specified  by  t;  therefore,  we  want  to  find  those 
values  of  t  which  correspond  to  the  t,  .   From  the  s, 
determined  by  interpolation  we  can  calculate,  using  the  law 
of  cosines,  a  set  of  corresponding  r,   values  and  from  them 
a  set  of  T]^  values.   We  can  then  find  the  values  of  S(T]^  ) 
by  interpolating  in  a  table  of  S(tj^),  where  t^  is  the  optical 
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depth  at  the  i  of  the  N  shells.  All  the  interpolations 
mentioned  above  are  performed  by  the  same  Aitkin-Lagrange 
interpolation  routine. 

The  table  of  N  values  of  S(tj^)  is  of  course  based  on 
the  previous  iteration  of  the  code. 

The  computation  to  determine  sets  of  T{r^)  ,    K(rj^), 
.S{ri,0,A),  I(rj^,9,A)  and  ultimately  F(rjg,A)  starts  from  a 
set  of  initial  values  of  T{r ■ ) ,  which  should  be  as  close  to 
the  expected  final  values  as  possible,  and  a  set  of  the 
presumed  physical  parameters.   The  set  of  K(r-j^)  can  be  found 
immediately  from  the  given  parameters  using  either  equation 
(5)  or  (7a).   Using  the  KCr^)  ,  the  initial  set  of  T(rj^)  and 
their  associated  optical  depths,  we  can  evaluate  the  integral 
in  equation  (23)  along  each  of  the  M  lines  radiating  from 
each  of  our  N  points.   Then  from  the  M  values  of  I(ri,9,A) 
obtained  about  each  of  the  N  points,  improved  sets  of  T  and 
S  can  be  calculated. 

This  process  is  repeated  until  successive  sets  of 
T(r^)  converge  to  within  a  few  tenths  of  a  degree.   The  sets 
of  T,  K,  S,  and  I  obtained  in  the  last  iteration  are  taken 
as  the  true  values  for  the  cloud.   A  set  of  F(rj^,A)  is 
obtained  from  the  I(r„,6,A)  using  equation  (8). 

To  save  computation  time  a  bolometric  case  with  u^  = 
0.0  is  run  to  determine  the  set  of  T(r^) .   If  in  the  numer- 
ical procedure,  we  are  performing  the  calculations  for  L 
values  of  A,  then  the  computation  time  will  be  cut  by  a 
factor  of  approximately  1/L  when  we  drop  the  A  dependence 
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and  run  a  bolometric  computation.   When  we  use  the  bolometric 
expressions  instead  of  their  wavelength  dependent  analogs, 
equation  (2  3)  becomes  the  bolometric  expression 

00 

I(r,e)^^^  =  /  e"'''  S[I(r,e)j^]  dt. 
o 

Like  equation  (23)  this  can  be  solved  iteratively  for  the 
specific  intensity.   The  set  of  T(rj^)  is  then  found  from 
equation  (22)  and  used  as  the  set  of  initial  values  in  the 
desired  wavelenth  dependent  solution. 

Note  that  we  can  assume  that  the  temperatures  deter- 
mined by  a  bolometric  computation  are  the  same  as  those 
that  would  result  from  a  detailed  many  wavelength  computa- 
tion only  if  the  particles  in  the  cloud  are  grey-bodies.   If 
we  want  to  consider  cases  where  Wq  =  cOq(A)  then  we  can  not 
use  a  bolometric  computation  to  determine  the  temperature 
distribution. 

The  basic  equations  in  Section  II  which  involve  inte- 
gration over  9  can  be  rewritten  such  that  the  integrations 
are  over  y  =  cos6  by  making  the  change  of  variable 
9  =  cos"  (u) ;  d9  =  -  dy/sin9  .    The  code  has  been  tested 
with  these  equations  programed  with  9  as  the  variable  of 
integration  and  with  cos9  as  the  variable  of  integration  - 
the  results  do  not  differ  significantly. 

When  the  central  star  can  not  be  considered  a  point 
source  it  is  necessary  to  integrate  over  the  angular  diameter 
of  the  star.   Since  the  specific  intensity  of  the  stars'  sur- 
face can  be  considered  constant  (we  can  ignore  the  small 
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effects  of  limb  darkening)  an  integration  from  9  =  0  to 
8  =  9*(rj^)  need  only  be  done  once  for  each  r^^  and  used  in 
place  of  the  relevant  point  -  source  approximation . 

For  a  shell  very  close  to  the  star  the  correct  contri- 
bution by  the  star  to  the  mean  intensity  and  the  approximate 
value  can  differ  by  as  much  as  a  factor  of  2.0  -  the  point  - 
source  approximation  giving  the  lower  value .   The  mean 
intensity,  J,  the  Eddington  flux,  H,  and  the  so-called  K- 
integral,  K,  are  all  moments  of  the  mean  intensity  I  (Mihalas, 

1970)  that  is,  M(n)  =  4  /■"■  I(u)y'^dvi, 

^  -1 
where  I  =  M(0) ,  H  =  M(l) ,  and  K  =  M{2) .   The  flux  F  should 

not  be  confused  with  H;  they  are  related  such  that  F  =  4TrH. 

When  RNR  >>  R*  we  need  not  worry  about  inaccuracies  due  to 

approximating  the  star  as  a  point  source.   For  example,  when 

R*  "^   0.2  X  lO-"-^  and  RNR  'v  0.15  x  10   we  find  an  error  due 

to  the  approximation  of  only  0.46%. 

Since  the  code  is  capable  of  evaluating  the  specific 

intensity  as  a  function  of  9  and  A  at  any  point  in  the  shell 

it  can  be  used  to  find  values  of  J,(t),  H  (t) ,  K, (x)  and  the 

Eddington  factor,  K^{t)/J^{t)  .   The  Eddington  factor  can 

vary  from  0.0  to  1.0  but  is  usually  assumed  to  be  1/3  (the 

Eddington  approximation,  K  =  1/3  J)  since  this  is  the  value 

it  would  have  in  an  isotropic  radiation  field. 


SECTION  IV 
RESULTS  OF  MODEL  CALCULATIONS 

The  computer  code  used  to  determine  the  radiation 
field  within  a  circumstellar  dust  shell  is  broken  into 
three  main  programs  (program  1,  program  2,  and  program  3). 
The  first  is  a  program  used  to  set  up  the  integration  grid 
and  store  it  in  the  computer  (usually  on  disk  storage) ,  the 
second  is  a  program  used  to  determine  the  bolometric  radia- 
tion field  within  a  dust  shell  and  to  find  the  temperature 
distribution,  and  the  third  program  is  used  to  calculate 
the  monochromatic  radiation  fields  for  a  number  of  wave- 
lengths using  the  temperature  distribution  determined  by 
the  second  program.   The  three  programs  must  initially  be 
run  in  the  order  given  but  program  1  need  only  be  run  once 
for  each  set  of  the  geometric  parameters  R^ ,  RNR,  R,  Tj, 
and  N.   Therefore  it  is  possible  to  run  a  number  of  differ- 
ent models  with  the  same  integration  grid  by  varying  the 
two  nongeometrical  parameters  T^  and  w^  . 

The  models  that  we  want  to  consider  here  are  charac- 
terized by  the  parameters  in  Table  1. 
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Table  1.   Values  of  the  seven  parameters  used  for  each  of 
the  models  computed  (cgs  units;  T*  in  °K) 


Model 

R* 

RNR 

R 

^I 

n 

'^^o 

T* 

I 

.202x10^3 

.15x10^^ 

.75x10^^ 

5.0 

1.5 

0.1 

6000 

II 

" 

II 

II 

0.0 

" 

III 

II 

II 

" 

1.0 

" 

IV 

" 

" 

". 

0.6 

II 

V 

" 

7.0 

" 

0.0 

" 

VI 

" 

II 

" 

0.1 

11 

VII 

II 

.202xl0-'-^ 

.202202x1013 

II 

" 

0.0 

" 

VIII 

" 

.625x10^5 

.75xl0-'-^ 

2.0 

2.0 

0.0 

II 

IX 

" 

.ISxlO^'* 

" 

5.0 

1.3 

0.1 

" 

X 

.8352xl0-'-^ 

II 

.8352x10^^ 

3.97 

0.0 

0.8 

12500 

One  additional  model,  model  XI,  will  be  described  later  in 
this  section. 

From  Table  1  we  can  see  that  the  most  common  values  of 
(i)  used  are  0.0  and  0.1.   This  is  because  the  amount  of 
computing  time  necessary  to  determine  the  radiation  field 
is  a  function  of  to  .   For  small  values  of  ai  ((jJq<.2)  program 
3  requires  approximately  six  iterations  to  determine  values 
of  the  flux,  for  example,  converged  to  three  significant 
figures,  whereas,  if  co^  =  1  the  program  will  need  at  least 
15  iterations  to  achieve  the  same  accuracy.   Because  of  the 
high  cost  of  models  with  large  values  of  cOq  only  two  are 
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computed,  models  IV  and  X.   Model  III  was  computed  for  only 
one  wavelength  and  therefore  is  not  comparable  in  cost  to 
models  IV  and  X  which  were  computed  for  11  wavelength.  When 
0)  =  1.0  we  can  save  time  by  running  the  program  for  only 
one  wavelength,  A-^say,  and  then  multiplying  the  quantities 
obtained  for  that  wavelength  by  the  factor  B(A2,T*)/B(A  ,T*) 
to  obtaine  the  quantities  appropriate  to  ^2 .   This  is  a 
valid  procedure  when  Wq  =  1.0  since  under  that  condition  all 
the  radiation  leaving  the  star  must  reach  the  surface  of  the 
shell  unchanged  in  its  spectral-energy  distribution. 

When  program  2  is  run  the  total  outgoing  bolometric 
flux  through  each  shell  is  computed  at  each  iteration.  Since 
the  only  source  of  radiation  is  the  central  star,  the  total 
outgoing  bolometric  flux  at  each  shell  should  be  equal  to 
the  bolometric  luminosity  of  the  star.   This  is  true  to  within 
0.5%  for  some  of  the  models  in  Table  1  and  to  within  4%  for 
all  the  models  in  Table  1  except  model  VII.   The  flux 
(bolometric)  emitted  by  the  outer  boundary  of  the  shell  in 
model  VII  is  too  low  by  23%  principally  because  of  the  very 
severe  geometry  used  in  model  VII.   This  model  was  run  as  a 
test  of  the  code  to  determine  its  accuracy  under  extreme 
conditions. 

To  determine  the  temperature  distribution  in  the  shell 
program  2  is  run.   The  temperature  distributions  for  the 
models  in  Table  1  are  shown  in  Figure  3.   Note  that  the 
curve  for  model  IX  is  very  nearly  a  straight  line.   This  is 
due  to  the  geometry  of  the  shell  and  the  fact  that,  as  we 
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shall  see,  the  radiation  field  in  the  shell  is  highly 
anisotropic.   For  a  highly  anisotropic  radiation  field  all 
moments  of  the  intensity  are  equal,  that  is  J  =  H  =  K.  Thus, 
J  «  F  and  we  know  that  T  <==  J-^^^   so  T  «  F  ''^  .   Because  of 
conservation  of  flux  in  a  spherical  shell  we  have  that 
F  cc  i/r2;  therefore,  T  ^c    (l/r^)^/'^  or  T  c^  l/r°'^.   There- 
fore, for  n  =  1.5  we  have  from  equation  (6)  that  T(r)  = 
C-,    +   C^/r^'^,  where  C^^  and  €3  are  constants.   This  means 
that 

T(l/r°-5)  =  Ci  +  C2  (l/r°-^) 
and,  therefore  we  can  see  that  t  will  be  a  linear  relation 
of  T. 

The  temperatures  of  the  dust  grains  as  mentioned  in 
Section  I  are  consistent  with  the  assumption  that  they  are 
carbon  grains  where  the  temperature  is  below  about  2100°K 
and  consistent  with  the  assumption  that  they  are  refractory 
silicates  for  temperatures  less  than  about  ISOO^K.   Note 
that  the  spread  in  temperatures  at  t  =  0  is  small  even 
though  some  of  the  models  have  very  different  physical 
parameters. 

Program  3  is  run  to  determine  a  variety  of  quantities 
for  the  cloud.   These  quantities  are:  I(X,T,e),  J(X,t), 
K(A,t)/J(X,t),  and  F(X,t).   Figures  4.  and  5  show  some 
representative  values  of  these  quantities  (computed  for 
model  IX) .   The  conversion  of  visual  to  infrared  radiation 
can  be  seen  on  Figure  4  by  following  the  relative  intensi- 
ties of  a  few  wavelengths  from  t  =  7  to  t  =  0 .   Figure  5 
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also  shows  this  by  way  of  the  graphs  for  J  and  F.   The  graphs 
for  K/J  show  that  the  radiation  field  is  generally  anisotropic 
although  at  depth  the  infrared  radiation  field  does  approach 
the  Eddengton  approximation,  K/J  =  1/3.   Thus,  the  use  by 
Huang  and  Apruzese  of  the  Eddington  approximation  can  be  ex- 
pected to  decrease  the  accuracy  of  their  results  significant- 
ly- 

In  Figure  6  the  run  of  temperature  with  optical  depth 
is  shown  for  two  models  calculated  with  this  code  and  for  two 
models  calculated  by  Huang  (19  71) .   The  two  models  calculated 
with  this  code  are  models  VIII  and  XI.   Model  XI  has  the  fol- 
lowing parameters:  x   =  2.32,   R*  =  .202  x  IQ-'"-^  cm, 
RNR  =  .609  X  lO"*-^  cm,   R  =  .75  x  lO"""^  cm,   n  =  2.0  and 
T*  =  6000 °K.   For  optical  depths  larger  than  about  2/3  the 
curves  for  models  VIII  and  XI  are  very  similar  the  curves  for 
Huang's  corresponding  models.    From  optical  depth  2/3  to  0 
the  curves  for  models  VIII  and  XI  drop  off  much  more  rapidly 
than  Huang's  curves.    This  is  just  what  we  would  expect  to 
happen  due  to  the  erroneous  value  of  the  Eddington  factor 
used  by  Huang  for  the  outer  part  of  the  cloud. 

Figure  7  shows  the  spectral-energy  distribution  for 
HD  45677  as  given  by  Swings  and  Allen  (1971)  and  Low  et  al. 
(19  70) .   On  the  same  graph  the  spectral-energy  distribution 
for  model  X  is  also  plotted.    Model  X  is  based  on  a  set  of 
parameters  derived  by  Apruzese  (1974)  for  HD  45677.    He 
determined  these  parameters  by  fitting  the  observations  in 
the  visual  since  his  model  is  not  capable  of  directly 
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fitting  the  infrared  observations  (actually  his  model  does 
use  the  infrared  data  in  that  he  attempts  to  make  the  total 
infrared  emission  calculated  by  his  model  equal  to  the 
total  infrared  emission  given  by  the  observations) .    The 
graphs  show  that  his  parameters  do  give  a  very  good  fit  in 
the  visual  region  and  that  the  total  emission  in  the 
infrared  predicted  by  his  parameters  is  approximately  cor- 
rect; however,  his  parameters  do  not  give  a  good  fit  to  the 
actual  infrared  spectral-energy  distribution.   A  different 
set  of  parameters  could  be  found  that  would  fit  both  the 
visual  and  the  infrared  regions. 

Figures  8  and  9  show  the  spectral  energy  distributions 
for  a  number  of  the  models  listed  in  Table  1.   On  Figure  8 
the  graphs  for  models  II,  V,  VIII  and  X  are  shown.   The 
curve  for  model  VI  would  fall  between  those  for  models  II 
and  V.   We  can  see  the  effect  of  changing  the  various  model 
parameters  by  studying  these  curves. 

The  curves  for  models  II  and  V  show  that  increasing 
the  total  optical  depth  of  the  shell  causes  the  peak  of  the 
distribution  to  move  toward  the  red  (Figure  3  shows  that  at 
the  same  time  the  temperature  of  the  inner  part  of  the  shell 
goes  up)  and  causes  the  amount  of  visual  radiation  filtering 
through  to  decrease  by  a  large  amount. 

Models  VIII  and  X  show,  as  we  would  expect,  that  when 
the  optical  depth  of  the  shell  is  decreased  by  a  large 
amount  the  central  star  causes  the  distribution  to  peak  in 
the  visual. 
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Figure  9  shows  the  spectral-energy  distributions  for 
models  I,  IV,  IX  and  the  means  of  the  observations  of  R  Mon 
due  to  Mendoza  (1966,  1968),  Low  and  Smith  (1966)  and  Low 
et  al.  (1970),  where  vertical  bars  show  the  amount  of 
variability  from  those  means.   The  vertical  error  bars  on 
the  observations  at  20  and  22 y  indicate  the  root-mean- 
square  error  of  a  single  observation,  approximately  20%. 

The  U  through  M  values  plotted  are  intensities 
calculated  from  Mendoza 's  observations  (corrected  for  inter- 
stellar extinction),  using  Johnson's  (1966)  calibration  of 
the  UBVRIJKLM  system.   The  interstellar  extinction  correc- 
tion was  made,  assuming  that  R  Monocerotis  is  located  in 
NGC  2264,  and  using  the  following  adopted  data  for  that 
cluster  and  hence  R  Monocerotis:  Eg_^  =0.06  (Strom  et  al. 
1971)  and  a  reddening  curve  of  the  NGC  22  44  type  with  R  =  5.4 
(Johnson,   196  8)  . 

Using  these  values  we  have  Ay  =  0.324,  and  from  the 
reddening  curve  we  obtain  the  extinction  in  magnitudes,  AM, 
for  each  of  the  observed  colors  (see  Table  2) . 

The  root-mean-square  errors  of  a  single  observation 
are  shown  for  the  20u  and  the  22y  observations  in  Figure  9, 
but  for  all  the  other  bands  the  root-mean-square  errors  of 
a  single  observation  are  negligible,  considering  the 
intrinsic  variability  of  the  object.   The  largest  such  er- 
ror is  for  the  M  band,  and  it  amounts  to  less  than  0.1 
magnitudes. 
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Table  2.   Interstellar  extinction  in  magnitudes.  Am,  for  R 
Monocerotis  in  each  of  the  colors  observed 


Band  Am 

U  0.424 

B  0.383 

V  0.324 

R  0.277 

I  0.224 

J  0.189 

K  0.141 

L  0.094 

M  0.083 

The  curve  for  model  IV  shows  that  a  large  increase 
in  the  value  of  cOq  causes  the  peak  of  the  distribution  to 
narrow  and  move  toward  the  blue,  it  also  shows  that  the 
direct  contribution  of  the  star  in  the  visual  increases 
greatly.   The  curves  for  models  I  and  IX  show  that  a  small 
increase  in  the  index  n  causes  the  peak  of  the  distribu- 
tion to  move  toward  the  blue,  while  the  ends  of  the  curve 
in  the  visual  and  far  infrared  drop. 

Model  I  is  the  best  fit  to  the  observations  of  R  Mon 
but  it  looks  as  though  Model  IX  with  a  slightly  larger  oj^ 
would  fit  somewhat  better. 

The  following  discussion  pertains  to  the  development 
of  the  model  for  R  Mon. 
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The  spectral  class  of  R  Monocerotis,  roughly  G  (Joy, 
1945),  and  Iximinosity  of  870  L^  (Low  and  Smith,  1966)  lead 
to  the  values  of  T*  and  R*  chosen  for  these  values  since 
4TTR*aT*^  ~  870L^.   The  value  for  x^  was  determined  from  a 
cursory  examination  of  the  observations  in  the  following 
manner.   We  have  the  ratio  of  the  visual  flux,  F^,  to  the 
bolometric  flux,  F_,  for  R  Monocerotis;  namely,  F^/Fg  -.006 

o 

(Low  and  Smith,  1966).   From  this  and  the  relation 
tj  =  iLn(FY/FB)  we  obtain  for  x  (Rj)  ;  Xj  ~  5. 

It  should  be  noted  that  this  method  gives  only  a 
lower  limit  to  xj,  since  we  cannot  say  what  portion  of  F^ 
is  actually  coherently  scattered  radiation.   Thus,  Wq  and 
T   cannot  be  determined  separately  with  any  ease.   The 
values  of  oo^,  Rj ,  R  and  n  used  were  just  rough  initial 
approximations . 

We  must  consider  the  following  observational  data 
for  R  Monocerotis:  (a)  the  spectrum  does  not  show  selective 
absorption  (Joy,  1945) ,  (b)  strong  variable  optical 
polarization  exists  indicating  scattering  by  large  grains 
(Zellner,  1970) ,  and  (c)  irregular  light  variations  occur 
of  several  magnitudes  (Joy,  1945;  Mendoza,  1966,  1968;  Low 
et  al.  1970).   In  addition  we  have  the  UBVRIJKLMQ  observa- 
tions of  Mendoza,  Low  and  Smith,  and  Low  et  al.  mentioned 
previously. 

Points  (a)  and  (b)  above  are  consistent  with  the  as- 
sumption that  the  cloud  is  made  up  of  solid  particles  which 
act  as  grey  absorbers  and  emitters.   For  a  period  of 
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variation  of  the  polarization  which  is  fairly  short,  the  as- 
sumption that  the  particles  are  randomly  oriented  should  be 
reasonable.   Because  of  the  variability  mentioned  in  point 
(c) ,  those  observations  in  several  colors,  which  were  made 
on  the  same  day  are  of  the  greatest  interest,  and  where  such 
observations  are  lacking,  we  are  forced  to  compare  the 
computed  fluxes  to  mean  observed  fluxes. 

Using  the  parameters  from  Case  I  given  by 

T^  =  5  , 

Rj  =  .15  X  10-^'^  cm  , 

R  =  .75  X  lO-'-^  cm  , 

n  =  1.5  , 

and  assuming  the  particles  to  have  the  density  of  carbon  or 

silicon  dioxide;  that  is,  p   ~  2.0,  we  have  from  equation 

(25)  for  the  mass  of  the  cloud. 

M  =  13dM   . 
® 

This  gives  us  for  particles  about  twenty  microns  in 
diameter  a  cloud  mass  of  about  10  ^K^    . 

The  set  of  parameters  for  Case  I  also  gives  us  for 
equation  (24)  the  following  form: 

9   ^_ 
p(r)  =  (1.5  x  10=")  ^1.5  • 

If  d  =  20u  and  r  ranges  from  lO"'-^  to  lO"""^  cm,  then  p(r) 
will  range  from  approximately  10"^  to  10"   grams/cm  . 
In  sximmary,  the  best  fit  to  the  observations  was 
obtained  by  a  model  with  the  following  parameters: 
Rj  =  .15  X  10^'*  cm  (~lAu)  , 
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R  =  .75  X  lO-"-^  cm  (-50  Au)  , 


T,  =  5  , 

n  =  1.5  , 
M  ~  10~^M 


lO"-"-^  <  P(r)  <  lO"-*-^  gm/cm^  , 


225°K  <  T  <  2,100°K  . 


Figure  1.   Radiation  is  shown  falling  on  point  1  at 
distance  r  from  the  central  star.   It  is 
incident  there  at  an  angle  8-,  after  being 
emitted  from  point  2  at  an  angle  9  .Part 
of  this  emitted  radiation  is  radiation 
from  point  3  which  was  incident  on  point 
2  at  an  angle  Qj'    ^^^   scattered  through 
an  angle  a.   Point  2  is  a  distance  r^ 
from  the  central  star  and  a  distance  s(t) 
from  point  1.   Analogous  triangles,  with 
vertices  at  the  star,  point  1  and  any 
point  on  the  line  directed  from  point  1, 
at  the  angle  0-,  ,  can  be  constructed. 

Note  that  the  angle  9,  is  designated 
simply  8  in  the  text. 
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STAR 


Figure  2.   A  schematic  diagram  of  the  integration  grid 
for  the  case  when  N  =  4 .  Four  lines  of  sight 
are  shown  to  indicate  the  manner  in  which 
the  boundary  conditions  at  r,  and  r^^  are 
handled.   The  filled  circles  on  the  lines 
represent  the  initial  grid  of  intersection 
points.   The  tick  marks  on  the  lines 
represent  the  points  found  at  equal  inter- 
vals of  T  by  an  iterative  technique.    The 
open  circles  represent  the  Gauss-Legendre 
points  for  a  three  point  quadrature  formula. 
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Figure  4.   Log  [I(A,e,T)]  vs.  9  for  model  IX.   The  angle 
9  is  in  radians.   The  three  panels  from  top 
to  bottom  are  for  t  =  7.0;  3.68;  0.0.   The 
curve  are  labeled  with  letters  designating 
their  respective  wavelengths  on  the 
UBVRIJKLMNQ  system.   Note  that  in  the  case 
for  T  =  7.0  the  star  is  not  considered  a 
point. 
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Figure  5.   Log  [J(A,t)],  K { A , t) /j( A , t) ,  and  Log  [F(A,t)] 
vs.  T  for  model  IX.   The  curves  are  labeled 
with  letters  designating  their  respective  wave- 
lengths on  the  UBVRIJKLMNQ  system. 
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